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PREFACE 


This  report  contains  a  user’s  manual  for  running  the  computer  program 
XYFREZ,  version  4.  The  basic  mathematical  method  incorporated  in  it  was  de¬ 
vised,  developed,  and  programmed  by  Dr.  Kevin  O’Neill,  of  the  Geotechnical  Re¬ 
search  Branch,  Experimental  Engineering  Division,  U.S.  Army  Cold  Regions  Re¬ 
search  and  Engineering  Laboratory.  The  initial  program  and  manual  were  ex¬ 
panded  and  tested  further  by  Dr.  O’Neill  and  Capt.  Richard  Jardine  of  the  U.S. 
Military  Academy  at  West  Point.  Support  for  this  work  came  primarily  from  the 
Department  of  Army  In-House  Laboratory  Independent  Research  Program  (ILIR). 
Documentation  and  work  to  make  the  program  more  “user  friendly”  was 
supported  in  part  by  DA  Project  4A762730AT42,  Task  BS,  Work  Unit  001 ,  Defor¬ 
mation  of  Soils  and  Pavements  Due  to  Freeze  !  Thaw ,  under  the  leadership  of  Dr. 
Richard  Berg. 

The  author  is  thoroughly  grateful  for  forward-looking  support  from  the  two 
sources  mentioned,  which  has  enabled  the  development  of  a  program  from 
innovative  mathematical  concepts,  all  the  way  through  to  a  documented  and 
reasonably  user-friendly  program.  The  diligent  and  insightful  assistance  of  Capt. 
Jardine  is  also  appreciated.  Terry  Rogers  and  Nancy  Humiston  reviewed  the 
report  with  great  care,  struggling  through  all  its  directions  in  detail,  as 
representative  users.  Their  reactions  were  much  appreciated  by  the  author,  and 
numerous  changes  were  made  in  response  to  their  experiences.  The  author  also 
acknowledges  the  patient  editorial  assistance  of  Steve  Bowen  and  Sandy  Smith, 
and  the  most  accommodating  illustrating  of  Ed  Perkins. 

The  contents  of  this  report  are  not  to  be  used  for  advertising  or  promotional  pur¬ 
poses.  Citation  of  brand  names  does  not  constitute  an  official  endorsement  or  ap¬ 
proval  of  the  use  of  such  commercial  products. 


ii 


CONTENTS 


Abstract  . 

Preface  . 

Introduction  . 

Input  data  specification 
Explanatory  comments 


Item  10  .  10 

Item  11  .  10 

Item  12  .  10 

Item  13  .  11 

Item  14  .  15 

Running  the  program  and  plotting  routines  .  16 

Sample  problem  .  17 

Blind  sample  problem  .  19 

Literatiire  cited  .  20 

Appendix  A:  Sample  problem  input  data  .  21 

Appendix  B:  Sample  problem  output  .  25 

Appendix  C:  Listing  of  the  program  .  37 


ILLUSTRATIONS 

Figure 

1 .  Exsunple  finite  element  mesh  constructed  of  linear  triangles  .  8 

2.  Example  boundary  condition  temperatures  over  time  .  13 

3.  Wedge-shaped  domain  of  the  sample  problem,  divided  into  linear 

triangles  as  plotted  by  KOPLOT  .  17 

4.  Cold  and  warm  side  boundary  temperatures  over  time  for  the 

sample  problem  .  18 

5.  Computed  points  for  freezing  firont  location  and  continuous  analyt¬ 

ical  solution  for  blind  seunple  problem  .  20 


Sion  For 

GFAtl 

TAB 

o.uiced 
.’Icpt,  ion_ 


Distribution/ 

1  Avallritillty  Codss 
j  Avail  Oiid/or 

iDlst  I  Special 


XYFREZ.4  User’s  Manual 


KEVIN  O’NEILL 


INTllODUCTION 

The  program  iises  linear  triangular  finite  elements  to  compute  two-dimensional 
heat  conduction  with  or  without  phase  change.  Convection  may  enter  through  con¬ 
vective  boundary  conditions,  as  explained  below,  but  is  not  considered  in  the  domain 
interior.  Heterogeneous  material  properties  may  be  specified.  The  current  version  is 
the  latest  in  a  developing  series,  each  with  various  added  efficiencies,  options,  and 
greater  idiot-proofing  than  the  previous.  The  user  should  note  that,  while  the  current 
version  has  been  tested  carefully  in  a  variety  of  problems  with  known  solutions,  it  has 
not  yet  been  widely  used.  This  causes  us  to  muse  on  an  old  adage — at  least  as  old  as 
the  yoimg  field  of  digital  computing — 

There  are  no  bug-free  programs,  only  programs  in  which  all  the  bugs 
have  not  yet  been  discovered. 

A  parallel  set  of  programs  is  being  developed,  the  RZFREZ  series,  for  solving 
axisymmetric  heat  conduction  and  phase  change  problems  in  cylindrical  coordinates. 

This  manual  describes  the  use  of  the  program  XYFREZ,  version  4.  Although  the 
program  utilizes  innovative,  state-of-the-art  methods  (CTNeill  1983)  both  the  program 
and  this  manual  have  been  constructed  as  much  as  possible  to  accommodate  users 
with  minimal  knowledge  of  computational  methods.  Ancillary  programs  for  auto¬ 
matic  generation  of  input  data  and  for  plotting  input  data  and  computed  results  are 
also  available  at  CRREL.  Contact  Kevin  (TNeill  with  any  questions  or  comments. 

The  mathematical  basis  of  the  methods  used  and  their  numerical  implementations 
are  described  in  the  paper  by  Kevin  O’Neill  (1983).  In  general,  an  explanation  of  the  fi¬ 
nite  element  method  and  its  use  in  the  program  is  beyond  the  scope  of  this  manual. 
However,  the  interested  reader  may  consult  any  of  a  variety  of  introductory  finite  ele¬ 
ment  books  to  learn  more  than  enough  to  grasp  this  program,  (e.g.  Finder  and  Gray 
1977,  and  Becker  et  al.  1981).  Sufhce  it  to  say  for  what  follows  that  one  approaches 
problem  solution  by  dividing  the  overall  physical  domain  up  into  a  mesh,  i.e.  smaller 
subsections,  each  of  which  is  known  as  a  finite  element.  In  the  particular  variant  of 


the  method  employed  here,  each  of  these  elements  must  be  a  triangle.  Based  on  linear 
interpolation  over  the  interior  of  each  element,  the  method  solves  for  temperature 
values  at  the  vertices  (“nodes”)  at  the  comer  of  each  element.  The  program  steps 
through  time  to  obtain  transient  solutions,  and  the  nodal  temperature  values  are 
printed  at  whatever  intervals  the  user  specifies.  The  number  of  elements  and  nodes 
that  are  used  in  a  given  problem  is  up  to  the  user,  who  inputs  the  mesh.  On  the 
whole,  more  elements  and/or  smaller  elements  can  be  used  where  one  wishes 
relatively  finer  resolution  to  be  applied  in  the  calculations.  An  example  mesh  appears 
in  the  section  entitled  Sample  Problem. 

In  principle,  you  may  use  as  many  nodes,  elements,  and  time  steps  as  you  please 
(or  can  afford)  in  a  given  problem.  The  program  is  set  up  so  that  specific  computer 
memory  requirements  are  derived  from  input  quantities,  while  a  large  block  of  space 
is  potentially  set  aside  in  the  introductory  routine  for  later  use  within  the  program.  If 
you  are  not  running  on  the  CRUEL  Prime  computer  system  or  do  not  have  a  virtual 
machine  at  your  disposal,  the  potential  memory  space  initially  set  aside  may  be  larger 
than  your  machine  can  handle.  For  example,  some  arrays  are  assigned  potential 
sizes  of  up  to  2000k  words.  If  you  cannot  live  with  that  kind  of  space  requirement  in 
“COMMON,”  caU  Kevin  O'NeiU. 

In  the  following  sections  the  required  input  data  are  described,  followed  by  ex¬ 
planatory  comments  on  the  meaning  of  and  options  within  that  data.  It  is  strongly 
recommended  that  you  study  the  explanatory  comments  carefully. 


INPUT  DATA  SPECIFICATION 

For  the  most  part,  data  for  this  program  are  read  in  firom  an  input  file  composed 
and  stored  by  the  user  before  running  the  program.  On  the  CRREL  Prime  system,  the 
program  will  ask  you  to  enter  the  name  of  yoxir  data  file  through  your  terminal,  after 
you  have  given  the  command  to  run  the  program.  On  other  systems  one  may  have  to 
alter  this  in  accordance  with  the  prevaihng  practices  and  available  options,  noting 
that  the  data  are  read  in  on  unit  5.  Unless  otherwise  noted,  everything  is  read  in  un¬ 
der  open  format.  That  is,  any  format  may  be  used  in  the  data  list — exponential  nota¬ 
tion  or  not,  commas  or  any  number  of  spaces  to  separate  entries,  etc.  The  usual  For¬ 
tran  variable  type  conventions  must  be  observed,  keyed  to  the  names  as  spelled  below: 
variable  names  beginning  with  the  letters  “P  through  “N”  are  integers  and  must  not 
contain  a  decimal  point.  Study  the  example  data  in  the  section  entitled  Sample 
Problem  for  more  guidance.  The  choice  of  unit  system  is  up  to  the  user,  but  must  be 
strictly  consistent. 


2 


, 


Item  No. 
0. 

Variable 

NBND,  NMBCl,  NELS, 

NTS,  MEDIA 

Meaning 

Number  of  nodes,  number  of  “first 
type”  boundary  nodes,  number  of 
elements,  number  of  time  steps, 
number  of  different  media  in 
region  simulated,  respectively 

1. 

10 

Unit  number  for  writing  all  output 
except  for  plotting  information 

2. 

TITLE 

Alphanumeric  output  title,  80 
chararacters 

3. 

ITRLM,  KPRTVL 

Iteration  limit,  printing  interval 
respectively 

4. 

IZDPLT,  MSHPRT, 

KBVPRT 

Integer  flags  for  printing  control 
of  output  for  isotherm  plotting, 
for  detailed  mesh  description,  and 
for  time  step  list  and  boundary 
condition  details,  respectively 

5. 

KPAUSE 

Flag  controlling  pause  option 

6. 

(NODE,  NDBC, 

XG,  YG) 

List  of  node  numbers,  boundary 
condition  indicator,  X  and  Y 
coordinates,  respectively 

7. 

(NL,  MTYPE,  Nl, 

N2,  N3) 

List  of  element  numbers,  material 
type  for  each  element,  and  the 
global  node  numbers  at  each 
element’s  vertices 

8. 

(HETF,  HETU,  CNDF, 

CNDU,  EL,  TF) 

List  of  thermal  properties  for 
each  MTYPE. 

9. 

TIMLIM,  AMP,  THETA 

Time  duration  of  simulation, 
numerical  control  parameters 

10. 

KIC 

Integer  flag  for  initial  condition 
choice  system 

11. 

TEMPIC  or 
(NODE,  TEMPIC) 

Initial  conditions  (see  explana¬ 
tory  comments) 

12. 

KBCT 

Choice  indicator  for  BC  and 
time  step  system  in  item  13 

13. 

(NODE,  TBCVAL),  and  KDT 
or 

(TIMBD,  NUMCHG,  NTSBD, 
(NODE,  TBCVAL)  ) 

List  of  boundary  nodes  and  their 
temperatures,  followed  by  time 
step  system  indicator 
(See  explanatory  comments) 

14. 


H,  TA 


Boundary  heat  transfer  coefficient, 
ambient  temperature 


EXPLANATORY  COMMENTS 


ItemO 

NBND  =  number  of  nodes  in  the  mesh,  including  the  boundary;  NELS  =  number  of 
elements  in  the  mesh;  NTS  =  the  number  of  time  steps  to  be  taken;  MEDIA  =  the 
number  of  different  media  (material  types)  within  the  region  covered  by  the  element 
mesh.  The  program  functions  by  solving  for  the  temperature  distribution  in  space  at 
sequentied  points  in  time.  It  proceeds  incrementally  from  each  point  in  time  to  the 
next,  hence  taking  “time  steps”  until  the  simulation  is  complete.  The  value  of  NTS 
only  indicates  the  total  number  of  such  steps  (points  in  time)  you  want  to  have.  The 
manner  in  which  those  steps  are  taken  is  determined  by  more  detailed  information  in 
other  items  in  the  input. 

NMBCl  =  the  number  of  boundary  nodes  with  prescribed  temperatures,  whether 
constant  or  time-varying  (as  opposed  to  boundary  nodes  with  flux  or  convective 
bovmdary  conditions).  These  are  “first  type”  boimdary  condition  nodes.  If  you  cannot 
easily  determine  this  number  or  in  any  case  want  to  have  the  program  do  so,  you  can 
begin  a  nm  using  any  integer  here.  The  program  determines  the  correct  value,  and  if 
your  input  is  wrong,  it  will  complain  mildly,  tell  you  the  correct  number,  and  then 
halt.  You  can  then  put  this  correct  number  in  your  data  and  start  over. 

Iteml 

The  user  supplies  an  integer  value  for  10  to  determine  where  the  output  will  be 
sent.  10  is  the  number  of  the  unit  on  which  all  non-plotting  output  will  be  printed.  On 
the  CRREL  Prime  system,  10  =  1  causes  the  output  to  appear  on  your  CRT  screen. 

Any  other  number  will  cause  writing  to  be  done  into  a  file.  If  you  specify  10  not  equal 
to  1,  the  program  will  ask  you  to  enter  an  output  file  name  through  your  terminal. 
This  output  file  name  may  correspond  to  a  new  file,  which  will  be  created,  or  an  old 
file  which  will  be  overwritten.  To  ensure  that  you  do  not  collide  with  Fortran 
preempted  file  numbers  or  others  defined  elsewhere  in  the  program,  use  10  =  6  to 
write  output  into  a  file. 

Users  on  systems  other  than  CRREL’s  Prime  may  have  to  rewrite  file  definitions, 
openings,  closings,  etc.,  in  accordance  with  the  conventions  built  into  their  0|}erating 
systems.  As  the  program  is  currently  set  up,  all  non-interactive  input  reading  is  done 
on  unit  5.  Unit  33  is  selected  for  writing  output  information  to  be  used  in  plotting  pro¬ 
grams  designed  for  the  CRREL  system. 


TITLE  is  an  80-character-long  label  for  your  problem.  In  other  words,  just  specify 
on  one  line  a  title  you  want  to  appear  at  the  top  of  the  output. 

Items 

ITRLM  is  the  iteration  limit.  Because  phase  change  problems  are  inherently  non¬ 
linear  and  because  we  wish  to  use  implicit  computational  schemes,  in  principle  some 
iteration  is  required  within  each  time  step  for  a  fully  consistent  solution.  Experience 
has  shown  however  that  this  is  usually  not  necessary  to  obtain  acceptable  answers, 
and  the  user  is  advised  to  specify  ITRLM  =  0,  at  least  for  starters.  If  time-step  sizes 
are  reasonable  (see  below)  but  bumpiness  still  develops  in  the  solution  in  either  space 
or  time,  higher  values  of  ITRLM  may  be  tried.  The  value  of  this  parameter  decides 
how  many  times  all  the  equations  will  be  set  up  with  improved  parameter  evaluations 
and  solved  over  again  (after  the  first  time)  within  a  single  time  step.  Thus  specifying 
higher  values  of  ITRLM  will  increase  computing  expenses,  assuming  that  time  step 
sizes  are  kept  the  same  as  ITRLM  is  increased.  Note  that  excessive  “bumpiness”  in 
the  shape  of  the  phase  change  isotherm  can  occur  when  it  is  curved  and  you  did  not 
specify  a  fine  enough  mesh  to  resolve  its  shape. 

KPRTVL  is  the  printing  interval.  That  is,  the  temperature  solution  for  all  nodes 
will  be  printed  out  every  KPRTVL  time  steps,  beginning  with  time  step  number 
KPRTVL. 

Item  4 

IZDPLT  is  an  integer  which  controls  printing  of  information  on  unit  33  for  later 
use  in  plotting.  As  of  this  writing,  the  plotting  program  to  do  this  has  only  been  imple¬ 
mented  on  the  CRREL  Prime  computer  system. 

If  IZDPLT  is  set  to  zero,  then  no  such  printing  is  done;  if  you  are  not  set  up  to  print 
on  unit  33  or  to  do  any  of  this  plotting,  set  IZDPLT  =  0. 

If  IZDPLT  =  -1,  then  only  your  input  data  for  the  mesh  will  be  stored  and  no  infor¬ 
mation  for  isotherm  plotting  will  be  saved.  This  option  allows  you  to  plot  your  input 
mesh  and  check  it  before  proceeding  with  actual  solution  of  the  problem  over  time.  See 
explanation  below  of  the  KPAUSE  parameter.  The  necessary  information  is  stored 
before  the  PAUSE  in  the  program  occurs.  Thus  you  can  halt  a  run  at  the  PAUSE 
point,  and  the  mesh  plotting  information  will  not  be  lost. 

If  IZDPLT  is  greater  than  zero,  then  information  will  be  saved  for  plotting  both  the 
mesh  and  fireezing  isotherms  on  it  at  specified  points  in  time.  Every  IZDPLT  time 
steps,  the  results  necessary  to  perform  that  plotting  are  dumped  into  a  file  you  have 
indicated,  beginning  with  time  step  number  IZDPLT. 


If  you  specify  any  nonzero  value  for  IZDPLT,  the  program  will  ask  you  to  enter 
through  your  terminal  a  file  name  under  which  plotting  information  for  the  current 
run  is  to  be  stored  for  later  use  (see  the  following  section).  The  name  may  be  chosen 
arbitrarily,  to  allow  you  to  distinguish  between  results  from  different  runs.  Later,  if 
you  use  the  program  KOPLOT  to  plot  the  contents  of  the  file,  then  that  program  will 
request  the  file  name. 

MSHPRT  is  a  flag  which  controls  printing  (on  unit  10)  of  mesh  information  which 
you  have  input.  If  you  specify  MSHPRT  =  1,  then  node  locations  together  with  NDBC 
values  will  be  written,  as  well  as  the  element  incidence  list  (which  nodes  go  with 
which  elements).  If  MSHPRT  =  0,  this  information  will  not  be  printed. 

The  value  of  KBVPRT  controls  the  printing  of  information  on  the  times  and  bound¬ 
ary  values  associated  with  each  time  step.  If  you  do  not  want  any  of  this  information 
printed,  specify  KBVPRT  =  0,  and  move  on. 

If  KBVPRT  =  1 ,  then  at  the  beginning  of  the  run  the  program  will  print  a  hst  of  the 
times  and  time  step  sizes  for  each  time  step,  without  indicating  the  associate  d  bound¬ 
ary  temperature  values. 

If  KBVPRT  =  2,  then  the  same  information  as  for  KBVPRT  =  1  will  be  printed,  to¬ 
gether  with  the  temperatures  of  all  t)rpe  1  boundary  nodes  at  each  time  step.  (See 
explanation  of  “type  1”  boundary  nodes  under  item  6  below.)  These  options  allow  you  to 
check  the  boundary  and  time  step  information  to  be  used  in  the  solution  procedures 
before  they  are  performed  (before  the  PAUSE).  This  is  especially  worth  checking  when 
you  have  had  the  program  interpolate  boundary  values  and  associated  times  from 
values  you  have  given  it  at  only  a  few  points  in  time  (see  items  12-13). 

Items 

This  item  concerns  the  PAUSE  option.  If  KPAUSE  =  1 ,  the  program  will  pause 
after  all  input  to  the  solution  procedures  has  been  either  input  by  the  user  or 
calculated  by  the  program;  you  will  then  be  asked  on  the  CRT  screen  (i.e.,  unit  1) 
whether  you  wish  to  proceed.  If  KPAUSE  =  0,  the  program  will  assume  that  you  do 
not  want  to  pause,  and  this  question  will  be  skipped.  This  option  is  intended  to  allow 
you  to  check  over  all  the  input  data  before  grinding  up  a  lot  of  computing  time 
marching  ahead  through  the  solution  in  time. 

Item  6 

Input  a  sequence  of  lines,  one  for  each  node,  in  which  the  node  number  (integer), 
the  boimdary  condition  indicator  NDBC  (integer),  and  the  X  and  Y  coordinates 
(floating  point)  for  each  node  are  given. 
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The  boundary  condition  indicator  NDBC  should  equal  1  if  the  associated  node  is  on 
the  boundary  and  temp)erature  values  are  to  be  specified  there  (“type  1”  boundary  con¬ 
ditions).  The  specified  temperature  may  be  constant  or  may  vary  with  time.  If  a  node 
is  on  the  boundary  and  NDBC  =  0  is  specified,  the  program  will  assume  that  the  node 
is  on  a  zero  flux  boundary  or  line  of  symmetry.  At  present  one  cannot  specify  nonzero 
boundary  flux  except  in  the  sense  that  applies  when  NDBC  equals  2,  explained  below, 
(jeneral  nonzero  flux  boundary  conditions  are  under  development. 

If  the  node  is  not  on  the  boundary,  specify  NDBC  =  0.  In  principle,  one  could  specify 
the  temperature  values  of  node  points  not  on  the  boundary,  as  explained  below  (items 
12-13)  by  simply  tagging  an  interior  node  with  a  value  of  1  for  NDBC.  Then  one  would 
provide  temperature  information  through  time  as  if  this  were  a  boundary  node.  The 
user  should  note  however  that  the  effect  of  constraining  interior  temperature  values 
has  not  been  investigated  in  the  development  of  this  program. 

If  NDBC  is  set  equal  to  2  for  a  node,  this  means  that  a  convective,  radiation,  or 
transfer-t3q)e  boundary  condition  applies  there.  This  means  that  heat  flux  throujth  the 
boundary  around  that  node  is  controlled  by  a  relation  of  the  form 

QN  =  H*{T-TA) 

where  QN  is  the  heat  flux  normal  to  the  boundary  (positive  outward  i,  H  is  n  (u  at 
transfer  coefficient  (heat  flow  rate/area/degree  ),  T  is  the  unknown  boundary  tt  niix  r 
ature,  and  TA  is  the  ambient  temperature.  Thus,  when  the  ambient  temperature  is 
lower  than  the  boundary  temperature,  heat  will  flow  out  of  the  region.  Values  for  // 
and  TA  are  given  by  the  user  in  item  14. 

Typical  input  for  item  6  might  look  like: 

1,  0,  0.0,  1.0 
8,  0,  10.0.  0.0 
2,  1,  3.0,  4.0 

These  example  numbers  say  that  node  number  2  has  a  “type  1”  boundary  condition  (its 
temperature  value  will  be  specified),  and  it  is  located  atX  -  3.0,  V  4.0.  Nodes 
number  1  and  8  will  not  have  their  temperatures  specified  and  thus  NDBC  =  0  for 
each.  They  are  either  interior  nodes  or  are  on  boundaries  where  the  normal  flux  of 
heat  through  the  boundary  (temperature  gradient  normal  to  the  boundary)  is  zero.  In 
either  case  the  program  will  solve  for  the  temperatures  at  nodes  1  and  8  and  they  are 
located  atX=  0.0,  T  =  1.0  andX  =  10.0,  T  =  0.0,  respectively. 


Figure  1.  Example  finite  element  mesh  constructed  of  lin¬ 
ear  triangles.  Shaded  elements  1  through  4  consist  of  ma¬ 
terial  number  2,  and  the  rest  are  material  number  1. 


Item? 

Specify  a  list  of  element  numbers,  with  each  one  followed  by  the  material  type  with¬ 
in  the  element,  in  turn  followed  by  the  node  numbers  of  each  vertex  of  the  triangular 
element.  For  example,  for  the  mesh  in  Figure  1  the  user  should  input; 


1,  2,  1,  4,  2 

2,  2,  2,  4,  5 

3,  2,  2,  5,  6 

4,  2,  2,  6,  3 

5,  1,  8,  4,  7 

6,  1,  5,  4,  8 

7,  1,  5,  8,  6 

8,  1,  8,  9,  6 


This  means  that  element  number  1  contains  material  number  2  and  has  nodes  1 ,  4, 
and  2  at  its  vertices,  while  element  number  5  contains  material  type  1  and  has  nodes 
8,  4,  and  7,  etc.  These  numbers  are  just  labels,  in  effect,  which  means  that  you  may 
number  the  elements  in  any  order  that  you  like.  Similarly,  you  can  number  the  nodes 
any  way  you  like.  The  relation  between  node  and  element  numbers  is  arbitrary,  as 
long  as  you  are  consistent  in  the  manner  shown  in  the  example  (i.e.,  make  sure  that 
you  associate  the  correct  nodes  with  each  element). 


The  efficiency  of  the  algorithm  for  ultimately  solving  the  equations  depends  greatly 
on  the  compactness  of  the  information  it  receives.  This  means  that  you  should  utilize 
your  freedom  in  node  numbering  to  minimize  the  difference  between  node  numbers 
across  any  one  element.  If  you  have  used  the  automatic  mesh  generator  program, 
then  your  node  numbering  scheme  will  already  have  been  optimized. 

If  you  have  the  program  print  out  detailed  mesh  information  (  MSHPRT  =  1),  you 
may  notice  that  the  nodes  associated  with  some  element  or  other  are  not  listed  in  the 
same  order  in  which  you  input  them.  That  is  because  the  program  automatically  re¬ 
orders  them  when  necessary  to  preserve  right-handedness  in  the  local  coordinate 
systems  in  which  it  operates.  (In  other  words,  nevermind.) 

Items 

This  item  consists  of  a  sequence  of  lines,  one  for  each  material  type  present,  which 
provide  thermal  properties  to  be  assigned  to  the  various  elements. 

Information  on  the  first  line  of  input  characterizes  material  number  1 ,  and  all  ele¬ 
ments  assigned  MTYPE  number  1  will  be  assigned  the  listed  characteristics.  The 
second  line  characterizes  material  number  2,  and  so  forth.  If  the  domain  considered 
is  entirely  homogeneous,  then  all  elements  would  have  been  assigned  MTYPE  =  1  in 
item  7,  and  item  8  would  consist  only  of  one  line. 

HETF  and  HETU  are  the  heat  capacities  of  the  frozen  and  unfrozen  material,  re¬ 
spectively,  per  unit  volume',  CNDF  and  CNDU  are  the  frozen  and  unfrozen  thermal 
conductivities,  respectively. 

EL  is  the  latent  heat,  meaning  the  latent  heat  released  per  unit  volume  of  material 
frozen.  Thus  for  pure  water  EL  would  be  the  latent  heat  per  unit  mass  of  water  times 
the  density  of  ice,  or  80  caFg  *  0.917  g/cc  =  73.36  cal/cc  ice.  For  soil,  this  figure  would  be 
multipUed  times  a  factor  equal  to  the  volume  of  ice  per  volume  of  soil.  For  a  rigid,  sat¬ 
urated  soil,  this  factor  equals  the  porosity. 

TF  is  the  phase  change  temperature. 

For  these  parameters,  as  for  everything  else  in  the  program,  any  consistent  system 
of  imits  may  be  used. 

Note  that  the  program  assigns  a  single  material  type  to  each  element,  but  allows 
for  a  discontinuity  in  material  (thermal)  properties  within  a  single  element  when 
phase  change  occurs. 

Item  9 

TIMLIM  is  the  length  of  time  to  be  simulated.  Its  value  is  ignored  smd  may  be 
specified  arbitrarily  when  boimdary  condition  information  is  given  for  specific  points 
in  time  (see  explanation  of  items  12  and  13). 


AMP  and  THETA  are  numerical  control  parameters  which  have  to  do  with  treat¬ 
ment  of  the  nonlinearity  and  the  degree  of  implicitness  in  the  solution  scheme,  re¬ 
spectively.  Recommended  values  are  0.65  for  each,  and  it’s  probably  best  not  to  mess 
with  these.  If  stability  problems  develop,  you  might  increase  THETA  to  0.75  or  to  1.0. 
AMP  is  an  interpolation  parameter,  which  determines  the  point  within  a  time  step  at 
which  nonlinear  coefficients  are  to  be  evaluated.  This  will  be  explained  in  a  forthcom¬ 
ing  paper.  In  general,  unless  you  know  what  you’re  doing,  leave  these  quantities 
alone.  The  value  of  AMP  only  matters  when  ITRLM  is  greater  than  zero. 

Item  10 

Here  the  program  is  told  bow  the  initial  conditions  are  going  to  be  specified. 

If  KIC  equals  zero,  then  all  nodes  have  the  same  initial  temperature,  to  be  specified 
in  item  11.  If  KIC  equals  1,  then  the  initial  temperature  for  each  individual  node  will 
be  specified  in  item  11. 

Item  11 

If  KIC  was  specified  as  zero,  then  a  single  vtdue  (TEMPIC)  is  input  here  and  all 
nodes  will  be  assigned  that  initial  temperature. 

If  KIC  was  specified  as  1 ,  then  this  item  contains  a  list  of  a//  the  nodes  with  the  ini¬ 
tial  temperature  of  each,  with  one  node  per  line.  For  example  the  list 

1,  10.0 
3.  4.0 
2,  9.55 

means  that  node  1  has  an  initial  temperature  of  10.0,  node  3  an  initial  temperature  of 
4.0,  and  node  2  an  initial  temperature  of  9.55. 

Item  12 

In  this  item  you  provide  the  program  with  a  value  for  KBCT,  which  indicates 
which  option  in  item  1 3  will  be  used  to  specify  the  boundary  conditions  and  to 
determine  the  time  steps.  KBCT  may  equal  either  1  or  2,  the  significance  of  each 
choice  being  explained  below  in  connection  with  item  13.  While  the  numerical  system 
used  here  has  been  designed  to  be  highly  stable,  it  is  always  possible  to  take  time  steps 
thoi.  are  so  large  that  they  cause  the  solution  to  degenerate.  A  conservative  rule  of 
thumb  is: 

Do  not  take  time  steps  that  cause  the  phase  change  isotherm  to  pass 
through  more  than  about  one  third  of  an  element  in  a  single  time  step. 


If  phase  change  is  not  occurring,  you  might  apply  this  to  some  other  reference  iso¬ 
therm,  an  appropriate  alternative  being  the  leading  edge  of  an  imposed  change.  You 
can  attempt  to  adhere  to  this  by  using  some  general  foreknowledge  of  the  problem  at 
hand  and  then  choosing  among  the  time  step  options  accordingly.  Alternatively  you 
can  try  some  crude  preliminary  runs,  and  subsequently  adjust  the  time  step  sizes. 

Item  13 

There  are  two  systems  of  boundary  condition  and  time  step  specification  under  this 
item,  each  of  which  pertains  only  to  setting  values  for  “t3rpe  1”  boundary  nodes.  That 
is,  we  are  only  talking  about  setting  temperatiire  values  for  certain  nodes,  in 
particular  those  for  which  NDBC  was  specified  as  1  in  item  6. 

BC,  DT  System  1 

If  KBCT  was  set  to  1 ,  this  means  that  all  type  1  boundary  nodes  will  have  tempera¬ 
tures  which  will  be  fixed  through  all  time.  Thus  item  13  contains  a  list,  with  one  node 
per  line,  of  all  t)q)e  1  nodes,  each  with  its  assigned  temperature. 

This  list  is  followed  by  a  value  for  KDT,  which  indicates  what  time  step  (DT)  system 
will  be  used.  If  KDT  is  assigned  the  value  1,  then  uniform  time  steps  are  taken  be¬ 
tween  time  zero  and  TIMLEM.  If  KDT  is  set  equal  to  2,  then  time  steps  are  taken 
which  are  equivalent  to  equal  increments  in  the  square  root  of  elapsed  time  between 
time  zero  and  TIMLIM.  The  reason  for  this  second  option  is  that  many  heat 
flow/phase  change  problems  with  constant  boundary  conditions  start  up  quickly  and 
then  settle  down,  with  changes  occurring  at  rates  more  or  less  proportional  to  the 
square  root  of  time.  In  any  case,  the  KDT  =  2  option  allows  you  to  take  small  time  steps 
in  early  time  with  increasingly  large  steps  as  time  proceeds. 

Note  again  the  danger  of  taking  time  steps  which  are  too  large,  which  can  creep  up 
on  you  unawares  if  you  just  let  the  program  project  DT  over  long  times.  The  implicit 
solution  system  used  is  designed  to  be  robust  in  the  face  of  large  time  steps,  even  when 
DT  is  so  large  that  accuracy  degenerates.  But  any  car  can  be  driven  off  the  road.  Ex¬ 
periments  with  the  program  have  also  suggested  that  the  numerical  method  used  is 
much  more  tolerant  of  large  time  steps  than  some  moving  boundary  approaches. 
Nevertheless,  the  nonlinearity  of  the  problem  undercuts  the  guarantees  of  stability 
derived  for  otherwise  comparable  implicit  linear  problems.  Investigations  of  moving 
boundary  phase  change  systems  suggest  that  excessively  laurge  time  steps  can  cause 
instability  even  when  there  is  very  little  action  in  the  solution.  The  most  common, 
seemingly  innocent  situation  where  one  is  likely  to  slip  into  unwise  DT  sizes  occurs 
near  a  steady  state.  As  time  runs  on,  little  happens  and  the  user  (or  the  program  on 
his  behalf)  is  tempted  to  take  larger  and  larger  time  steps.  If  instability  should  develop 
in  the  solution,  check  this  item. 
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The  KDT  time  step  options  are  summarized  as  follows: 

KDT  =  1  Uniform  time  steps 

KDT  =  2  Uniform  steps  in  square  root  of  time 

If  there  were  only  two  boundary  nodes  (numbered  3  and  1 2),  both  at  8  degrees 
through  all  time,  and  you  wanted  to  take  uniform  time  steps,  you  would  specify  the 
following  for  item  13,  assuming  you  had  specified  KBCT  =  1  in  item  12: 

3,  8.0 
12,  8.0 
1 

BC,  DT  System  2 

If  KBCT  =  2,  you  must  specify  values  for  each  of  the  type  1  boundary  nodes  at  vari¬ 
ous  times.  When  you  specify  the  time  at  which  some  boundary  values  apply,  the  pro¬ 
gram  will  also  record  how  many  time  steps  you  want  to  be  taken  since  the  last  such 
specification.  Thus  in  addition  to  boundary  values  you  will  be  indirectly  specifying 
time  step  sizes  and  overall  time  duration.  This  option  for  boundary  condition  and  time 
step  specification  may  seem  complicated  in  its  verbal  description.  However,  with  a 
little  familiarity  you  will  see  that  it  is  designed  to  give  you  the  greatest  flexibility  in 
assigning  time  and  space  dependent  boundary  temperatures  with  arbitrary  time 
steps,  by  specifying  a  minimum  of  information.  The  system  is  most  easily 
comprehended  through  an  example,  which  will  be  followed  by  a  more  general 
statement  of  the  system. 

Suppose  you  had  only  two  type  1  boundary  nodes,  numbered  1  and  4,  and  you 
wanted  to  assign  their  temperatures  and  take  time  steps  as  shown  in  Figure  2.  In¬ 
formation  need  only  be  specified  for  the  points  indicated  by  the  arrows  on  the  horizon¬ 
tal  axis  in  the  figure.  The  values  at  time  zero  should  not  be  specified  here — they  have 
already  been  specified  via  the  initial  conditions.  For  the  total  elapsed  time  TIMLIM  of 
30  you  would  input  the  information  in  Figure  2  as 

10.0,  2,  3 
1,  5.0 

4,  1.5 
30.0,  1,  2 
1,  5.0 

This  means  that  the  first  boundary  information  is  to  be  specified  for  an  elapsed  time  of 
(TIMED  =)  10.0.  Two  nodal  temperatures  will  be  specified  (NUMCHG  =  2),  and 
(NTSBD  =)  three  time  steps  will  be  taken  since  the  last  information  was  given  (in  this 
case,  since  the  initial  condition  shown  on  the  left  axis  of  the  figure).  The  next  two 
lines  contain  the  specific  boundary  temperature  information  for  time  10.0:  the  values 
5.0 


Figure  2.  Example  boundary  condition  temperatures  over  time.  Vertical  arrows  indi¬ 
cate  points  in  time  where  information  is  given  in  the  input  data. 


and  1.5  are  assigned  to  nodes  1  and  4,  respectively.  Because  three  time  steps  have 
been  indicated,  the  program  will  take  three  time  steps  after  time  zero,  each  of  equal 
length  (10.0/3  =  3.33  =  DT)  and  all  boundary  temperatures  will  be  linearly  interpolated 
over  the  intermediate  time  steps  since  the  initial  condition.  This  is  indicated  on 
Figure  2  between  time  zero  and  the  first  vertical  arrow. 

The  next  line  gives  information  for  an  elapsed  time  of  30.0.  Only  the  temperature  of 
one  node  is  to  be  assigned,  and  2  time  steps  are  indicated  since  the  last  specification. 
Boundary  values  are  given  in  the  next  line  for  time  30.0,  the  temperature  value  of 
node  number  1  is  specified  as  5.0.  The  program  then  determines  uniform  time  steps 
between  time  10.0  and  time  30.0,  and  linearly  interpolates  the  temperature  values  over 
intermediate  time  steps. 


\*  v'A.’ 


1^'^  ti*  fi*  ii*  ti*  A*  ft*  S*  V  Vi^i*  Vdti*  ia*  li‘  l^  r^  f>‘  I 


Note  that  for  the  second  point  in  time  only  the  temperature  of  node  1  had  to  be 
specified.  Temperature  information  need  only  be  given  when  the  rate  of  change  of  a 
boundary  temperature  is  to  be  altered.  Because  the  rate  of  change  of  node  4  remained 
unaltered  after  time  10.0,  nothing  was  indicated  for  it  at  time  30.0,  and  the  program 
extrapolated  along  the  same  line  begun  over  the  first  three  time  steps.  In  effect,  a 
linear  interpolation  is  accomplished  between  time  zero  and  time  30.0.  This  system 
allows  you  to  hold  the  temperature  of  a  boundary  node  constant  by  assigning  it  the 
same  value  in  two  sequential  specifications,  and  then  giving  no  more  information 
about  it.  So,  for  example,  if  you  wanted  to  continue  holding  node  1  at  5.0  degrees 
during  time  steps  after  time  step  5  (after  time  30.0),  you  could  simply  omit  any  infor¬ 
mation  about  it  in  subsequent  specifications  after  those  shown  above. 

The  temperature  of  a  type  1  boundary  node  will  be  kept  at  its  initial  condition  until 
and  unless  some  information  to  the  contrary  appears  in  this  section  of  the  input.  For 
example,  suppose  there  were  a  type  1  boundary  node  which  you  wanted  to  remain  at 
its  initial  temperature  while  nodes  1  and  4  changed  as  above.  In  this  case  you  would 
input  exactly  the  same  information  for  times  10.0  and  30.0  as  appears  above.  If  you 
wanted  the  temperature  of  this  additional  boundary  node  to  begin  changing  at  time 
30.0,  you  would  specify  a  new  value  for  it  for  some  time  after  30.0. 

At  first  it  may  seem  strange  that  boundary  temperatures  must  be  indicated  when  a 
nodal  temperature’s  rate  of  change  is  altered,  and  not  when  its  value  is  altered.  This 
is  to  allow  efficient  assignment  of  boundary  temperature  changes  which  are  smooth 
(linear)  over  relatively  long  stretches  (e.g.  that  of  node  4  in  the  example)  without  the 
user  having  to  indicate  every  single  change  in  value.  If  the  program  were  set  up  to  al¬ 
low  specification  of  values  by  node  (i.e.  not  by  time)  using  only  beginning  and  ending 
values  over  arbitrary  time  intervals,  that  would  require  the  user  to  coordinate  the  in¬ 
formation  in  a  potentially  complex  way,  so  that  the  same  time  steps  fit  the  changes  in 
boundary  values  for  all  nodes. 

One  always  retains  the  option  of  inputting  values  for  all  type  1  boundary  nodes  at  all 
time  steps  (or  every  so  many  time  steps).  This  is  in  fact  likely  to  be  suitable  for  labora¬ 
tory  or  field  data,  in  which  data  are  often  obtained  in  (long)  columns  of  parallel  mea¬ 
surements  over  the  entire  duration  of  am  operation.  While  it  might  be  possible  to 
achieve  the  saune  results  by  specifying  less  information,  there  is  no  compelling  com- 
putationail  reason  for  doing  so. 

If  you  only  wamt  to  change  the  time  step  size  but  let  the  program  continue  with  the 
previously  specified  boundauy  vailues,  simply  specify  NUMCHG  =  0  at  a  new  time,  and 
provide  no  further  information.  That  is,  in  the  example  above,  if  after  time  30.0  you 
wamted  node  1  to  continue  at  5  degrees  and  node  4  to  continue  changing  at  the 
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same  rate,  but  you  wanted  to  take  four  more  time  steps  of  DT  =  25.0,  you  would  follow 
the  specifications  above  with  this  one: 

130.0,  0,  4 

The  overall  system  when  KBCT  =  2  can  be  stated  briefly  as  follows: 

Boundary  temperature  values  for  at  least  some  type  1  nodes  (NDBC  =  1 )  must  be 
read  in  for  some  time  or  times  after  the  initial  condition.  Thus  when  KBCT  =  2,  item 
13  consists  of  a  sequence  of  boundary  information  specifications.  Each  specification 
consists  of  a  value  for  TIMBD  (the  time  for  which  information  is  being  given),  a  value 
for  NUMCHG  (the  number  of  nodes  whose  rates  of  change  or  initial  values  are  to  be 
altered),  and  a  value  for  NTSBD  (the  number  of  time  steps  to  be  taken  since  the  last 
specification).  This  information  is  followed  as  part  of  the  same  specification  by  a  list 
NUMCHG-long  of  nodes  each  followed  by  its  temperature  at  TIMBD.  The 
temperature  of  a  type  1  boundary  node  will  remain  at  its  initial  condition  through  all 
points  in  time  for  which  no  other  information  is  provided. 

Note  that  no  value  ofKDT  is  to  be  provided  when  KBCT  =  2. 


Item  14 

In  this  item  one  specifies  values  for  a  heat  transfer  coefficient  and  ambient 
temperature  for  those  boundary  nodes  flagged  with  NDBC  =  2.  If  NDBC  does  not  equal 
2  for  any  node,  then  skip  this  item. 

At  present,  the  program  is  set  up  so  that  only  uniform  constant  values  for  H  and 
for  TA  apply  for  all  such  boundary  nodes.  This  should  be  generalized  shortly  so  that 
variable  H  and  TA  can  be  used,  and  so  that  H  can  be  a  function  of  the  temperature. 

A  potential  problem  arises  at  nodes  where  different  types  of  boundary  conditions 
intersect.  To  keep  the  rules  straight  here,  one  might  think  of  a  hierarchy  of  boundary 
conditions,  with  the  highest  being  a  type  1  condition  (NDBC  =  1 ),  the  second  highest  a 
heat  transfer/radiative/convective  type  (NDBC  =  2),  and  the  lowest  being  a  specified 
(zero)  flux  (NDBC  =  0).  The  higher  priority  BC  type  always  takes  precedence.  For 
example,  with  reference  to  Figure  1,  suppose  that  you  wanted  the  temperature  to  be 
specified  along  the  boundary  line  containing  nodes  1,  2,  and  3,  and  you  wanted  a 
convective-type  condition  over  tbe  line  between  nodes  3,  6,  and  9.  In  this  case  node  3 
would  be  a  type  1  boundary  node,  and  you  would  specify  its  temperature  eilong  with 
that  of  nodes  2  and  1 .  Alternatively,  if  the  side  1 ,  2,  3  was  a  zero-flux  boundary,  and  the 
side  2,  6,  9  a  convective  boundary,  then  node  3  would  be  treated  as  a  convective¬ 
boundary  node.  If  one  observes  these  conventions,  then  the  program  is  always  capable 
of  translating  boundary-node  information  correctly  into  information  pertaining  to 
boundary  areas  between  the  nodes. 
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There  is  currently  one  constraint  on  the  use  of  convective-type  boundary  conditions: 
Do  not  specify  all  three  nodes  of  a  single  element  as  being  convective-boundary  nodes. 
Thus,  if  one  wanted  a  convective  condition  to  apply  over  the  boundary  between  nodes  4, 
7,  and  8,  in  Figure  1  he  would  have  to  use  a  different  mesh  arrangement.  The  pro¬ 
gram  is  set  up  to  alert  you  if  you  have  violated  this  rule. 

RUNNING  THE  PROGRAM  AND  PUOTTING  ROUTINES 

If  you  are  not  computing  on  the  CRREL  Prime  system,  you  will  have  to  compile  the 
source  code,  and  make  sure  in  the  process  that  all  the  conventions  of  your  system  are 
satisfied,  including  the  availability  of  input/output  files.  The  remaining  paragraphs 
in  this  section  refer  only  to  operations  on  the  CRREL  system. 

When  you  want  to  run  XYFREZ  on  the  CRREL  Prime  system,  first  make  sure  that 
you  have  stored  all  necessary  input  in  a  file  in  your  current  user  file  directory.  The 
file  may  have  any  name.  To  execute  the  program  simply  enter  the  command 

XYFREZ 

as  if  it  were  a  Prime  system  command.  The  program  will  begin  executing  as  soon  as 
that  command  has  registered,  and  you  will  promptly  be  asked  the  name  of  your  data 
file.  Enter  that  name  and  the  program  will  proceed,  writing  information  either  on  the 
screen  or  in  a  file,  depending  upon  what  you  have  specified  for  10. 

During  program  execution,  a  plotting  output  file  may  be  generated,  depending  up¬ 
on  the  value  you  have  specified  for  IZDPLT.  At  CRREL,  this  file  will  reside  in  your 
user  file  directory  until  you  remove  it.  You  can  use  the  plotting  software  described 
below  to  obtain  mesh  and  phase  change  isotherm  plots  for  the  run  which  generated 
the  file  in  question. 

To  plot  the  mesh  only  (for  IZDPLT  =  -1  in  your  input)  or  to  plot  both  the  mesh  and  a 
number  of  phase  change  isotherm  locations  (IZDPLT  greater  than  zero),  you  may  run 
the  program  KOPLOT.  The  program  resides  on  the  CRREL  system  and  is  designed 
for  it.  It  may  also  work  for  comparable  systems,  or  for  individual  office  setups,  but 
that  must  be  determined  on  a  case-by-case  basis. 

To  nm  the  program  KOPLOT  on  the  CRREL  Prime  system,  simply  assign  the  HP 
plotter  in  the  terminal  room  to  yoxu^lf  by  entering  the  command 

ASHP 

After  the  system  has  reassured  you  (“OK”),  then  enter 

KOPLOT 


When  that  command  has  registered,  the  plotting  program  will  begin.  It  will  ask  you 
to  enter  the  name  of  the  file  in  which  your  plotting  data  are  stored,  i.e.,  you  will  have 
to  enter  the  name  you  supplied  in  response  to  the  question  asked  by  XYFREZ. 
KOPLOT  is  currently  set  up  so  that  isotherms  are  labelled  by  time  step,  rather  than  by 
time,  and  the  graphs  are  8^2  by  11  inches  in  size.  An  example  is  provided  in  the  next 
section. 


SAMPLE  PROBLEM 

This  section  presents  a  sample  problem  to  demonstrate  the  input,  execution, 
and  graphics  associated  with  XYFREZ.  The  example  is  a  simple  case  for  which  there 
exists  a  known  solution.  The  numerical  solution  of  the  problem  using  XYFREZ  was 
originally  presented  in  the  reference  by  OTJeill  (1983).  A  listing  of  the  input  data  is 
provided  in  Appendix  A. 
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Figure  3.  Wedge-shaped  domain  of  the  sample  problem,  divided  into  linear  triangles 
as  plotted  by  KOPLOT.  Dashed  lines  show  freezing  front  locations,  with  associated 
time  step  numbers. 
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Figure  4.  Cold  and  warm  side  boundary  temperatures  over  time  for  the  sample  • 

problem. 


Graphical  description  of  the  boundary  conditions,  mesh  geometry,  and  the  location 
of  the  phase  transition  in  time  are  contained  in  Figures  3  and  4.  Note  that  the  domain 
in  question  is  the  shape  of  a  wedge,  in  which  the  straight  boundaries  are  formed  by 
radial  lines  of  symmetry.  The  normal  flux  of  heat  is  zero  along  those  boundaries,  and 
NDBC  is  set  to  zero  there.  The  nodes  on  the  left  and  right  ends  of  the  wedge  have 
NDBC  =  1 ,  and  their  temperatures  will  be  specified  as  part  of  the  boimdary  condition 
input  as  per  Figure  4.  The  domain  is  described  in  terms  of  two  material  types,  so  that 
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one  may  see  the  input  operations  which  that  entails.  Both  materials  are  given  the 
same  characteristics,  however,  so  that  this  corresponds  to  the  original  problem  for 


homogeneous  material,  as  reported  by  (TNeill  (1983). 


Note  in  the  data  that  bouikbuT  taB^wratures  are  giyen  below  for  nodes  1-3  and  22- 
24  every  other  time  step.  Hie  program  intorpolatee  the  times,  time  steps,  and  bound¬ 
ary  values  in  between.  The  ou^mt  listed  below  can  be  made  more  compact  if  one 
specifies  zero  for  MSHPRT  and  KBVPRT.  The  computed  results  in  Appendix  B 
below  correspond  to  those  reported  by  (yNeOl  (1983)  and,  as  may  be  seen  in  that  paper, 
are  quite  accurate.  Freezing  firout  Ineatjona  are  shown  on  the  mesh  in  Figure  3,  as 
plotted  by  the  KOPLOT  routine,  the  nee  of  iHuch  is  described  above. 

It  is  recommended  that  the  user  eiqwriment  with  XYFREZ  using  this  sample  data 
before  proceeding  to  a  real  problem.  If  you  have  obtained  this  program  firom 
Kevin  O’Neill,  you  probably  have  this  izqmt  data  on  tape  or  in  a  computer  file.  Tty 
altering  input  temperatures  or  time  steps,  diange  the  material  properties  of  patches 
of  elements,  etc.,  and  note  the  effects  of  doing  so. 

BUND  SAMPLE  PROBLEM 

Information  is  included  in  this  section  f<^  a  "blind”  sample  problem,  to  help  the 
user  gain  experience  in  using  this  program.  That  is,  data  are  given  in  somewhat  the 
same  form  in  which  they  might  be  available  for  a  real  problem,  in  terms  of  physical 
parameter  values,  length  of  time  to  he  simulated,  boundary  conditions,  etc.  However, 
the  user  must  decide  on  and  specify  the  particulars  of  mesh,  time  step  size,  and  other 
matters  of  judgment  or  individual  strategy  required  to  run  the  program.  The  problem 
poeed  is  a  simple  one,  with  an  exact  solution,  which  is  also  provided. 

The  problem  consists  of  one-dimensional  breezing  in  a  semi-infinite  medium.  The 
initial  temperature  is  4.0  degrees,  the  cold-side  boimdary  temperature  is  constant  at 
-10.0  degrees,  and  the  freezing  temperature  is  0.0  degrees.  Volumetric  heat  capacities 
in  the  firozen  and  unfirozen  zones  are  0.49  and  0.62,  respectively.  Frozen  and  unfrozen 
thermal  conductivities  are  9.6E-3  and  6.9E-3  respectively.  Hie  latent  heat  content  is 
17.68  per  imit  voliime  firozen  material.  These  values  correspond  roughly  to  parame¬ 
ters  which  might  characterize  a  moist  soil,  in  cal-cgs-deg  C.  The  exact  solution  to  this 
problem  is  given  by 

Xf=:B*V’t 

where  t  is  time,  XfiB  the  location  of  the  breezing  firont  relative  to  the  cold  side,  and  B  is 
a  constant  equal  to  8.8865E-2. 
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Figure  S.  Confuted  points  for  freezing  front  location  and  continuous  analytical 
solution  for  blind  sample  problem. 

Figure  5  shows  niunerical  locations  of  the  freezing  firont  versus  the  continuous 
analytical  solution,  over  40  time  steps,  up  to  a  time  of  1 .0E!4.  Of  course,  your  own 
answers  will  differ  somewhat,  depending  on  your  choices  for  time  step,  mesh 
characteristics,  and  ITKLM  value.  You  will  have  to  obtain  freezing  firont  locations  by 
interpolating  between  nodes  on  either  side  of  the  firont,  based  on  their  temperatures. 
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APPENDIX  A:  SAMPLE  PROBLEM  INPUT  DATA 
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o.oo(«-ai 
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3.827E+01 
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1 
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18 
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19 

0 
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20 

0 
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22 

1 
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23 
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1 
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4 

1 
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5 

1 
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7 

1 
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8 
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9 
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10 

1 
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11 

1 

9  8  11 
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1 
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13 

1 
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0.4900EX»  6  2 

1  -3.240E^0 
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0.5625B+09  6  2 

1  -3.393E+00 
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22  1.543E^ 
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24  1.413EUOO 
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XYFREZ  Sample  Problem:  Freezing  of  a  Wedge 
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28 
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BC 
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4.000E+00 

2 

4.000E+00 
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6 

4.000E+00 
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8 
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9 
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10 
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11 

4.000E+00 
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13 
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14 

4.0O0E-f00 

15 
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16 
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18 
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their  temperatures 

1  3.250E+00 

2  3.250E+00 

3 

3.250E+00 

22  4.000E+00 

23  4.000E+00 

24 

4.000E+00 

TIME  STEP  =  2 

Nodes  and 
1  2.501E+00 

22  4.000E+00 


TIME  =  2.500E+06 

their  temperatures 
2  2.501E+00 

23  4.000E+00 


DT  =  1.250E+06 

3  2.501E+00 

24  4.000E+00 


TIME  STEP  =  3  TIME  =  6.250E+06 


DT  =  3.750E+06 


Nodes  and  their  temperatures 
1  1.780E+00  2  1.780E+00 

22  3.999E+00  23  3.999E+00 


3  1.780E+00 

24  3.999E+00 


TIME  STEP  =  4 

Nodes  and 
1  1.060E+00 

22  3.998E+00 


TIME  =  l.OOOE+07 

their  temperatures 
2  1.060E+00 

23  3.998E+00 


DT  =  3.750E+06 

3  1.060E+00 

24  3.998E+00 


TIME  STEP  =  5 


TIME  =  1.625E+07 


DT  =  6.250E+06 


Nodes  and 
6.102E-01 


their  temperatures 
2  6.102E-01 


1 


3  6.102E-01 


22  3.971E+00 


TIME  STEP  = 


23  3.971E+00 


TIME  =  2.250E+07 


Nodes  and  their  temperatures 
1.605E-01  2  1.605E-01 

3.944E+00  23  3.944E+00 


TIME  STEP  = 


3.125E+07 


Nodes  and  their  temperatures 


-1.575E-01 

3.870E+00 


TIME  STEP 


2  -1.575E-01 

23  3.870E+00 


TIME  = 


4.000E+07 


Nodes  and  their  temperatures 
1  -4.754E-01  2  -4.754E-04 

22  3.796E+00  23  3.796E+00 


TIME  STEP  = 


TIME  =  5.125E+07 


Nodes  and  their  temperatures 


-7.212E-01 

3.692E+00 


2  -7.212E-01 

23  3.692E+00 


TIME  STEP 


6.250E+07 


Nodes  and  their  temperatures 


1  -9.669E-01 

22  3.588E+00 


TIME  STEP  *  11 


-9.669E-01 

3.588E+00 


7.625E+07 


Nodes  and  their  temperatures 
1  -1.168E+00  2  -1.168E+00 


24  3.971E+00 


DT  =  6.250E+06 


3  1.605E-01 

24  3.944E+00 


DT  =  8.750E+06 


3  -1.575E-01 

24  3.870E+00 


DT  »  8.750E+06 


3  -4.754E-01 

24  3.796E+00 


DT  -  1.125Et-07 


3  -7.212E-01 

24  3.692E+00 


DT  -  1.125E+07 


3  -9.669E-01 

24  3.588E-«-00 


DT  «  1.375E+07 


3  -1.168E+00 


m  ^ 


22  3.473E+00 


TIME  STEP  =  12 


23  3.473E+00 


TIME  =  9.000E+07 


Nodes  and  their  temperatures 
1  -1.369E+00  2  -1.369E+00 

22  3.359E-t-00  23  3.359E-t-00 


TIME  STEP  =  13 


TIME  = 


1.062E+08 


Nodes  and  their  temperatures 
■1.539E+00  2  -1.539E+00 

3.245E4-00  23  3.245E+00 


TIME  STEP  =  14 


TIME  =  1.225E+08 


Nodes  and  their  temperatures 
1  -1.709E+00  2  -1.709E+00 

22  3.130E-t-00  23  3.130E  +  00 


TIME  STEP  =  15 


TIME  =  1.412E+08 


Nodes  and  their  temperatures 
1  -1.857E+00  2  -1.857E+00 

22  3.019E+00  23  3.019E+00 


TIME  STEP  *  16 


TIME  = 


1 .600E+08 


Nodes  and  their  temperatures 
1  -2-004E+00  2  -2.004E+00 

22  2.908E-t-00  23  2.908E+00 


TIME  STEP  =  17 


1 .812E+08 


Nodes  and  their  temperatures 
1  -2.134E-t-00  2  -2.134E+00 


24  3.473E+00 


DT  =  1.375E+07 


3  -1.369E+00 

24  3.359E+00 


DT  =  1.625E+07 


3  -1.539E+00 

24  3.245E+00 


DT  =  1.625E+07 


3  -1.709E+00 

24  3.130E+00 


DT  =  1.875E+07 


3  -1.857E+00 

24  3.019E+00 


DT  =  1.875E+07 


3  -2.004E+00 

24  2.908E-t-00 


DT  =  2.125E+07 


3  -2.134E+00 


2.803E+00 


23 


2.803E+00 


24 


2,803E+00 


TIME  STEP  »  18 

TIME  *  2.025E+08 

DT 

-  2.125E+07 

Nodes 

-2.264E+00 

and 

their  temperatures 

2  -2.264E+00 

3 

-2.264E+00 

2.698E+00 

23  2.698E+00 

24 

2.698E-t-00 

TIME  STEP  =  19 

TIME  =  2.262E+08 

DT 

=  2.375E+07 

Nodes 

-2.381E+00 

and 

their  temperatures 

2  -2.381E+00 

3 

-2.381E+00 

2.599E+00 

23  2.599E-t-00 

24 

2.599E+00 

TIME  STEP  =  20 

TIME  =  2.500E+08 

DT 

-  2.375E-t-07 

Nodes 

-2.497E+00 

and 

their  temperatures 

2  -2.497E+00 

3 

-2.497E+00 

2.501E+00 

23  2.501E+00 

24 

2.501E+00 

TIME  STEP  =  21 

TIME  =  2.762E+08 

DT 

=  2.625E+07 

Nodes 

-2.602E+00 

and 

their  temperatures 

2  -2.602E+00 

3 

-2.602E+00 

2.408E-«-00 

23  2.408E+00 

24 

2.408E+00 

TIME  STEP  -  22 

TIME  =  3.025E+08 

DT 

*  2.625E+07 

Nodes 

-2.707E+00 

and 

their  temperatures 

2  -2.707E+00 

3 

-2.707E+00 

2.316E+00 

23  2.316E+00 

24 

2.316E+00 

TIME  STEP  -  23 

TIME  -  3.312E+08 

DT 

»  2.875E+07 

Nudes 

-2.803E-f00 

and 

their  temperatures 

2  -2.803E+00 

3 

-2.B03E+00 

^  ■»- 


22  2.229E+00 


TIME  STEP 


23  2.229B4^00 


TIME 


3.600t«0S 


Nod«a  and  thair  taapara  tnrM 
-2.899E+00  2  >2.899B4>00 

2.142E+00  23  2.142B‘*>00 


TIME  STEP 


TIME 


3.912B<H)8 


Nodes  and  their  teaperatures 
1  -2.988E-*-00  2  >2.98884-00 

22  2.061E+00  23  2.061B4-00 


TIME  STEP 


4.225E4-08 


Nodes  and  their  temperatures 
1  -3.076E-J-00  2  -3.076E-400 

22  1.979E+00  23  1.979E+00 


TIME  STEP  -  27 


TIME 


4.562E-408 


Nodes  and  their  temperatures 
1  -3.158E+00  2  -3.158E-t-00 

22  1.902E+00  23  1.902E-*-00 


TIME  STEP  »  28 


TIME  »  4.900E4-08 


Nodes  and  their  temperatures 
1  -3.240E+00  2  -3.240E4-00 

22  1.825E-»-00  23  1.82SE400 


TIME  STEP 


TIME  -  5.262B408 


Nodes  and  their  teaperatures 
1  -3.317E-t-00  2  -3.3I7E4-00 


v -V-'W 


24  2.2298^ 


Vt  •  2.8758-407 


3  -2.89984^ 

24  2.14284-00 


DI  -  3-125E-407 


3  -2.98884-00 

24  2.061E400 


DT  -  3.125E-407 


3  -3.07684-00 

24  1.97984-00 


DT  -  3.375E-K)7 


3  -3.15884-00 

24  1.90  284-00 


DT  -  3.3758-407 


3  -3.2408-400 

24  1.82584-00 


DT  -  3.6258-407 


3  -3.317E+00 


l 


22  1.753E+00 


TIME  STEP  =  30 


23  1.753E+00 


TIME 


5.625E+08 


Nodes  and  their  temperatures 
1  -3.393E+00  2  -3.393E+00 

22  1.680E-»-00  23  1.680E-«-00 


TIME  STEP  =  31 


6.012E+08 


Nodes  and  their  temperatures 
1  -3.464E+00  2  -3.464E+00 

22  1.612E+00  23  1.612E+00 


TIME  STEP  =  32 


TIME 


6.400E-t-08 


Nodes  and  their  temperatures 
1  -3.535E+00  2  -3.535E+00 

22  1.543E+00  23  1.543E+00 


TIME  STEP  =  33 


TIME  =  6.812E+08 


Nodes  and  their  temperatures 
1  -3.602E+00  2  -3.602E+00 

22  1.478E+00  23  1.478E+00 


TIME  STEP 


TIME 


7.225E+08 


Nodes  and  their  temperatures 
1  -3.669E+00  2  -3.669E+00 

22  1.413E+00  23  1.413E+00 


TIME  STEP 


TIME 


7.662E+08 


Nodes  and  their  temperatures 
1  -3.732E+00  2  -3.732E+00 


24  1.753E+00 


DT  =  3.625E+07 


3  -3.393E+00 

24  1.680E+00 


DT  »  3.875E+07 


3  -3.464E-t-00 

24  1.612E+00 


DT  =  3.875E+07 


3  -3.535E+00 

24  1.543E+00 


DT  =  4.125E+07 


3  -3.602E+00 

24  1.478E+00 


DT  =  4.125E+07 


3  -3.669E+00 

24  1.413E+00 


DT  »  4.375E+07 


3  -3.732E+00 


32 


1.352E+00 

23  1.352E+00 

24  1.352E+00 

TIME  STEP  =  36 

TIME  -  8.100E+08 

DT  =  4.375E+07 

Nodes 

-3.795E+00 

1.290E+00 

and 

their  tenperatures 

2  -3.795E+00 

23  1.290E+00 

3  -3.795E+00 

24  1.290E+00 

TIME  STEP  -  37 

TIME  >  8.562E-t-08 

DT  =  4.625E+07 

Nodes 

-3.855E+00 

1.231E+00 

and 

their  tenperatures 

2  -3.8S5E-»-00 

23  1.231E+00 

3  -3.855E+00 

24  1.231E+00 

TIME  STEP  =  38 

TIME  -  9.025E-»-08 

DT  =«  4.625E+07 

Nodes 

-3.915E+00 

1.172E+00 

and 

their  teaperatures 

2  -3.915E+00 

23  1.172E+00 

3  -3.915E+00 

24  1.172E+00 

TIME  STEP  =  39 

TINE  >  9.512E-»-08 

DT  -  4.875E+07 

Nodes 

-3.971E+00 

1. 116E+00 

and 

their  temperatures 

2  -3.971E+00 

23  1.116E-»-00 

3  -3.971E+00 

24  1.116E-t-00 

TIME  STEP  =  40 

TIME  -  l.OOOE+09 

DT  -  4.875E+07 

Nodes 

-4.028E+00 

1 .060E+00 

and 

their  temperatures 

2  -4.028E+00 

23  1.060E+00 

3  -4.028E+00 

24  1.060E+00 

Calculated  Teaperatura  Solution 


Plotting  data  atorad  for  TIHB  *  1.625B-t-07 


TIME  - 

6.250B-t>07 

TIME  STEP 

=  10 

NODE 

TEMP 

NODE 

TEMP 

NODE 

TEMP 

1 

-9.669B-01 

2 

-9.669B-01 

3 

-9.669E-01 

4 

-8.441B-02 

5 

-8.567E-02 

6  ' 

-8.364E-02 

7 

5.637E-01 

8 

5.616E-01 

9 

5.634E-01 

10 

1.450B-»-00 

11 

1.448E+00 

12 

1.450E+00 

13 

2.058E-»-00 

14 

2.053E+00 

15 

2.058E+00 

16 

2.681B-«-00 

17 

2.674E+00 

18 

2.681E+00 

19 

3.207E+00 

20 

3.200E+00 

21 

3.207E+00 

22 

3.588E+00 

23 

3.588E+00 

24 

3.588E+00 

Plotting 

data  stored 

for 

TIME  = 

6.250E+07 

Plotting 

data  stored 

for 

TIME  - 

1.412E+08 

TIME  = 

2. 

500E+08 

TIME  STEP  * 

20 

NODE 

TEMP 

NODE 

TEMP 

NODE 

TEMP 

1 

-2.497E+00 

2 

~2.497E+00 

3 

-2.497E+00 

4 

-1.612E+00 

5 

-1.612E+00 

6 

-1 .611E+00 

7 

-9.806E-01 

8 

-9.813E-01 

9 

-9.809E-01 

10 

-1.017E-01 

11 

-1.025E-01 

12 

-1 .021E-01 

13 

5.515E-01 

14 

5.500E-01 

15 

5.511E-01 

16 

1.259E+00 

17 

1.257E+00 

18 

1.259E+00 

19 

1 .924E-«-00 

20 

1.921E+00 

21 

1 .924E+00 

22 

2.501E+00 

23 

2.501E-*-00 

24 

2.501Et-00 

Plotting 

da  ta 

stored 

for 

TIME  - 

2.500E+08 

Plotting 

da  ta 

stored 

for 

TIME  - 

3.912E+08 

TIME  -  5.625E+08  TIME  STEP  -  30 


TEMP 

NODE 

TEMP 

NODE 

TEMP 

-3.393E+00 

2 

-3.393E+00 

3 

-3.393E+00 

-2.498E+00 

5 

-2.498E+00 

6 

-2.497E+00 

-1 .859E+00 

8 

-1.859E+00 

9 

-1 .860E+00 

-9.658E-01 

11 

-9.658E-01 

12 

-9.662E-01 

-3.284E-01 

14 

-3.289E-01 

15 

-3.287E-01 

3.907E-01 

17 

3.894E-C1 

18 

3.906E-01 

1 •068E+00 

20 

1.067E+U0 

21 

1 .OeSE+OO 

1 .680E+00 

23 

1.680E+00 

24 

1.680E+00 

Plotting 

data  stored 

for  TIME  = 

5.625E+08 

Plotting 

data  stored 

for  TIME  = 

7.662E+08 

TIME  = 

l.OOOE+09 

TIME  STEP  = 

40 

TEMP 
-4.028E+00 
-3. 134E+00 
-2.495E+00 
-1.601E+00 
-9.636E-01 
-2.594E-01 
4.356E-01 
1.060E+00 


NODE  TEMP 

2  -4.028E+00 

5  -3.134E+00 

8  -2.495E+00 

11  -1.601E+00 

14  -9.638E-01 

17  -2.615E-01 

20  4.357E-01 

23  1.060E+00 


NODE  TEMP 

3  -4.028E+00 

6  -3.133E+00 

9  -2.495E+00 

12  -1.602E+00 

15  -9.639E-01 

18  -2.594E-01 

21  4.355E-01 

24  1.060E+00 


Plotting  data  stored 


for  TIME  = 


l.OOOE+09 


36 


APPENDIX  C:  LISTING  OF  THE  PROGRAM 


00001 : 
00002: 
00003: 
00004: 
00005: 
00006: 
00007: 
00008: 
00009: 
00010: 
00011: 
00012: 
00013: 
00014: 
00015: 
00016: 
00017: 
00018: 
00019: 
00020: 
00021 : 
00022: 
00023: 
00024: 
00025: 
00026: 
00027: 
00028: 
00029: 
00030; 
00031 : 
00032: 
00033: 
00034: 
00035: 
00036: 
00037 ; 
00038: 
00039; 
00C40: 
00041 : 
00042: 
00043 : 
00044 : 
00045: 
00046: 
00047 ; 
00048: 
00049: 
00050 : 
00051 : 
00052 : 
00053 : 
00054: 


PROGRAM  XYFREZ.4.PUBLIC3 


This  program  solves  the  2-D  heat  equation,  including  phase  change 
effects  through  the  heat  capacity. 

Main  Program 


This  routine  obtains  input  for  determining  array  dimensions. 
The  "real"  main  program  follows  in  SUBROUTINE  MAIN 

COMMON  NDBC(2000),  XG(2000),  YG(2000), 

&  R(2000),  T(2000),  TOLD(2000),  TLAST ( 2000 ), NGLNBC ( 2000 ) , 

&  ANS(2000),  MTyPE(4000),  NBCELG ( 1000 ) ,  NBCNGL ( 1 000 ) , 

&  N0DEL( 20000 ) , DT( 2000) ,HETU( 10) , HETF( 10) , 

&  CNDU{ 10 ) , CNDF( 10) , TF( 10) / EL( 10) ,  A ( 200000 ), TBC ( 2000000 ) 


CHARACTER*32  DATFIL 


WRITE( 1 , 701 ) 

701  FORMAT(///////,5X, '  WELCOME  TO  XYFREZ,  VERSION  4',/////, 
&  'Enter  name  of  data  file') 

READ(1,7777)  DATFIL 
7777  FORMAT(A) 

OPEN( 5 , FILE=DATFIL) 


READ(5,*)  NBND,  NBCIIN,  NELS,  NTS,  MEDIA 


IF(NBND.LE. 2000)  GO  TO  110 
WRITE( 1 , 7001 )  NBND 

7001  F0RMAT( /////,' NBND  value  of ',110, 'is  greater  than  2000!') 
GO  TO  1000 

no  IF(NBC1IN,LE.1000)  GO  TO  120 
WRITE( 1 , 7002 )  NBCIIN 

7002  FORMAT! /////,' NMBCl  value  of, 110, 'is  greater  than  1000!' 
GO  TO  1000 

120  IF  I NELS . LE. 4000 )  GO  TO  130 
WRITE! 1  ,  7003  )  NELS 

7003  FORMAT! /////,' NELS  value  of, 110,  'is  greater  than  4000!') 
GO  TO  1000 

130  IF ! NTS . LE . 2000 )  GO  TO  140 
WRITE! 1 , 7004 )  NTS 

7004  FORMAT! /////,' NTS  value  of', 110, 'is  greater  than  2000!') 
GO  TO  1000 

140  IF(MEDIA.LE. 10 )  GO  TO  150 
WRITE! 1  ,  7005  )  MEDIA 

7005  FORMAT! /////,' MEDIA  value  of', 110, 'is  greater  than  10!') 


m 


00055 

00056 

00057 

00058 

00059 

00060 

00061 

00062 

00063 

00064 

00065 

00066 

00067 

00068 

00069 

00070 

00071 

00072 

00073 

00074 

00075 

00076 

00077 

00078 

00079 

00080 

00081 

00082 

00083 

00084 

00085 

00086 

00087 

00088 

00089 

00090 

00091 

00092 

00093 

00094 

00095 

00096 

00097 

00093 

00099 

00100 

00101 

00102 

nni03 

00104 

00105 

00106 

00107 

00108 


1000 

7999 


WRITE( 1 , 7999 ) 

FORMAT(/////////, '  EXECUTION  HALTED’ 
CLOSE(5) 

STOP 


)  CONTINUE 

NTSP  =  NTS  +  1 

The  following  is  a  dummy  value  of  NEW/  for  use  in  the 
dimension  statement  in  SUBROUTINE  MAIN 

NEW  =  1 

CALL  MAIN(A,  NBND,  NEW,  NDBC,  XG ,  YG ,  R,  T,  TOLD,TLAST, 
&  NGLNBC,  ANS,  MTYPE,  NBCFLG,  NBCIIN,  NBCNGL,  DT, 

&  NTSP,  TEC,  NODEL,  NELS,  MEDIA) 


SUBROUTINE  MAIN(A,  NBND,  NEW,  NDBC,  XG ,  YG ,  R,  T,  TOLD, 

&  TLAST, NGLNBC ,  ANS,  MTYPE,  NBCFLG,  NBCIIN,  NBCNGL,  DT ,  NTSP, 
&  TBC,  NODEL,  NELS,  MEDIA) 


DIMENSION  A(NBND,NBW) ,  NDBC(NBND),  XG(NBND),  YG(NBND) 

DIMENSION  R(NBND),  T(NBND),  TOLD(NBND),  TLAST(NBND) 

DIMENSION  NGLNBC(NBND) ,  ANS(NBND),  MTYPE(NELS),  NBCFLG ( NBCl IN ) 
DIMENSION  NBCNGL ( NBCl IN ) ,DT( NTSP) , TBC (NTSP, NBCl IN) , NODEL ( NELS , 3 ) 
DIMENSION  TITLE(20),  TF(IO),  EL(IO) 

DIMENSION  HETUdO  )  ,  HETF(IO)  ,CNDU(  10)  ,  CNDF(  10) 


CHARACTER*32  OUTFIL,  PLTFIL 


NTS  =  NTSP  -1 


INPUT 

READ(5,*)  10 
IF(IO.EQ.l)  GO  TO  12 
WRITE( 1 ,773) 

FORMAT(///,'  Enter  file  name  for  general  output',///) 
READ(1,7777)  OUTFIL 
FORMAT( A) 

OPEN(IO,FILE=OUTFIL) 

CONTINUE 


>  '  «  'J  V  'J' 


00109:  READ(5,1791)  ( TITLE ( K ), K= 1 , 20 ) 

00110:  1791  FORMAT(20A4) 

00111:  WRITE( 10, 791 )  ( T ITLE ( K ) , K= 1 , 20 ) 

00112:  791  FORMAT(/////20A4//////////////) 

00113:0 

00114:  READ(5,  *)  ITRLM,  KPRTVL 

00115:  WRITE(IO,707)  NBND,  NBCIIN,  NELS,  NTS,  ITRLM,  KPRTVL 

00116:  707  FORMAT(/,10X,  'NBND'  ,5X,  'NMBCl'  ,6X,  'NELS'  ,7X,  'NTS'  , 

00117:  &  5X, ' ITRLM' ,4X, 'KPRTVL* ,/,4X, 6110) 

00118:0 

00119:  READ(5,*)  I ZDPLT , MSHPRT , KBVPRT 

00120:  WRITE(IO,719)  I ZDPLT , MSHPRT , KBVPRT 

00121:  719  FORMAT!//, 21X, ' IZDPLT' ,6X, 'MSHPRT' ,6X, 

00122:  &  'KBVPRT'  ,/,15X, 4112) 

00123:0 

00124:  IF( IZDPLT.EQ.O)  GO  TO  112 

00125:  WRITE(1,774) 

00126:  774  FORMAT!////,'  Enter  file  name  for  plotting  output',///) 

00127:  READ!1,7777)  PLTFIL 

00128:  0PEN!33,FILE=PLTFIL) 

00129:  112  OONTINUE 

00130:0 

00131;  WRITE! 1 ,8001 ) 

00132:8001  FORMAT !///////, 20X , '  Preparing  Finite  Element  Mesh...',/////) 

00133:  READ!5,*)  KPAUSE 

00134:  WRITE! 10, 709)  KPAUSE 

00135:  709  FORMAT! //, 30X, ' KPAUSE  =’,I2) 

00136:0 

00137:0 

00138:0  The  next  loop  reads  each  node,  the  type  be  at  the  node 
00139:0  and  the  coordinates  of  the  node.  Also  sets  be  flags. 

00140:0 

00141:  NMBOl  =  0 

00142:  KB02  =  0 

00143:  DO  100  NODE  =  1 , NBND 

00144:  READ!5,*)  ND , NDBOV , XGV , YGV 

00145:  IF!NDBGV.EQ.2 )  KBG2  =  1 

00146:  NGLNBG!ND)  =  0 

00147:  IF!NDBGV.NE. 1 )  GO  TO  95 

00148:  NMBOl  =  NMBOl  +  1 

00149:  NGLNB0!ND)  =  NMBCl 

00150:  NBCNGL! NMBCl )  =  ND 

00151:  95  NDBC!ND)  =  NDBOV 

00152:  XG!ND)  =  XGV 

00153:  YG!nd)  =  YGV 

00154:  100  CONTINUE 

00155:0 

00156:  IF!NMBC1 .EQ. NBCIIN)  GO  TO  111 

00157:  WRITE!1,7770)  NBCIIN,  NMBCl 

00158:  7770  FORMAT! ///,  lOX ,' NMBCl  value  of, 16,'  input  is  different  from' 
00159:  &  /,10X,'  actual  value  of, 110, 

00160:  &  ///////, 20X, 'EXECUTION  HALTED!') 

00161:  CLOSE!5) 

00162:  CLOSE! 10) 
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00163: 
00164: 
00165: 
00166:0 
00167:0 
00168: 
00169: 
00170: 
00171 : 
00172  :  0 
00173:0 
00174:0 
00175:0 
00176: 
00177  : 
00178: 
00179: 
00180: 
00181 : 
00182: 
00183:0 
00184:0 
00185:0 
00186:0 
00187:0 
00188:0 
00189 
00190 
00191 
00192 
00193 
00194 
00195 
00196 
00197 
00198 
00199 
00200 
00201 
00202 
00203 
00204 
00205 
00206 
00207 
00208 
00209 
00210 
00211 
00212 
00213 
00214 
00215 
00216 


CLOSE( 33  ) 
STOP 

111  CONTINUE 


IF(MSHPRT.EQ.l  ) 

&  WRITE( 10, 701 )  (K,NDBC(K) ,XG(K) , YGiK) ,K=1 ,NBNU 
701  FORMAT{//, (12X, 'Node' ,4X, •  BC  ',12X,'X', 

&  lOX, ' Y' ) ,/, (9X,2I7,3X,1P2E11.2)  ) 

The  following  loop  reads  the  connectivity  of  the  nodes 
and  determines  the  material  type  of  the  element 

DO  110  NLMT  =  1,NELS 

READ(5,»)  NL,LTYPE,N1,N2,N3 
MTYPE(NLMT)  =  LTYPE 
NODEL(NL,l)  =  N1 
N0DEL(NL,2)  =  N2 
NODEL(NL,3)  =  N3 
110  CONTINUE 


Check  for  flipped  element:  is  X1*Y2  -  Y1*X2  .LT.  0 
Also  determine  NDG  and  NBW 
NNDG  =  0 

DO  250  L  =  1,NELS 
KNTERR  =  0 

249  N1  =  N0DEL(L,1) 

N2  =  N0DEL(L,2) 

N3  =  N0DEL(L,3) 

NLBC2  =  NDBC(Nl)  +  NDBC(N2)  +  NDBC(N3) 

IF(NLBC2.NE.6)  GO  TO  2249 
WRITE( 10, 7249  )  L 

7249  FORMAT( //// , '  Element  no', 15,'  has  three  type  2  BCs',// 

&  'FIX  IT! ' ) 

GO  TO  1111 

2249  NDIF  =  IABS{N2-N1) 

IF(NDIF.GT.NNDG)  NNDG  =  NDIF 
NDIF  =  IABS(N3-N2) 

IF{NDIF.GT.NNDG)  NNDG  =  NDIF 
NDIF  =  IABS(N1-N3) 

IF{ NDIF. GT. NNDG)  NNDG  =  NDIF 
XI  =  XG(N2)  -  XG(Nl) 

X2  =  XG(N3)  -  XG(Nl) 

Y1  =  YG{N2)  -  YG(Nl) 

Y2  =  yG{N3)  -  YG(Nl) 

PROD  =  X1*Y2  -  Y1*X2 
IF(PROD.GT.O.O)  GO  TO  250 
KNTERR  =  KNTERR  +  1 
IF{KNTERR.GT. 1  )  THEN 
WRITE! 10,7111)  L 

FORMAT!////,'  Incurable  mesh  problem  in  element  no'. 


7111 


GO  TO  1111 


00217: 

00218:  ELSE 

00219:  ENDIF 

00220:  NHOLD  =  NODEL(L,3) 

00221:  NODEL(L,3)  =  NODEL(L,2) 

00222:  NODEL(L,2)  =  NHOLD 

00223:  GO  TO  249 

00224:  250  CONTINUE 

00225:C 

00226:C  Arguments  to  be  passed  for  later  calculations  involving  the  matrix  k, 
00227:C  e.g.  NEW  is  the  bandwidth  of  the  banded  matrix. 

00228:C 

00229:  NEW  =  2*NNDG  +  1 

00230:  NDG  =  NNDG  +  1 

00231  :C 

00232:C 

00233:C 

00234:  IF(MSHPRT.NE. 1  )  GO  TO  1250 

00235:  WRITE( 10,702)  ( K , MTYPE ( K ) , ( NODEL { K , J ) , J=1 , 3 ) , K=1 , NELS ) 

00236:  702  FORMAT (//, 2 ( 6X Element 3X Ma ter ia 1 6X , '  Nodes 2X ),/ , 

00237:  &  2(9X, I4,5X, I2,6X,3I4)  ) 

00238:C 

00239:  1250  CONTINUE 
00240:C 

00241:C  Write  mesh  plotting  info  on  disk 
00242:C 

00243:  IF(IZDPLT.EQ.O)  GO  TO  1108 

00244:  WRITE(33,7010)  NELS 

00245:  7010  FORMAT(I7) 

00246:  DO  107  NL  =  1,NELS 

00247:  WRITE(33,7011)  NL,  MTYPE(NL),  ( NODEL( NL , J ) , J=1 , 3 ) 

00248:  7011  FORMAT(5I7) 

00249:  107  CONTINUE 

00250:  WRITE(33,7010)  NEND 

00251:  DO  108  ND  =  1 , NEND 

00252:  WRITE(33,7012)  ND,  NDEC(ND),  XG(ND),  YG(ND) 

00253:  7012  FORMAT ( 2 17 , 1P2E14 . 5 ) 

00254:  108  CONTINUE 

00255:  1108  CONTINUE 

00256:C 

00257:C  Input  thermal  properties 
00258:C 

00259:  DO  1109  1=1, MEDIA 

00260:  READ(5,*)  HETF ( I ) , HETU ( I ) , CNDF { I ) , CNDU ( I ) , EL ( I ) , TF ( I ) 

00261 :  WRITE( 10,705)  I, HETF ( I ) ,HETU( I ) ,CNDF( I ) ,CNDU( I ) , EL ( I ) , TF ( I ) 

00262:  705  FORMAT( ///, 1 5X ,' Thermal  Properties  of  Material 13 ,//, 9X , 

00263:  &  'HETF' ,6X, 'HETU' ,6X, 'CNDF' ,6X, ’CNDU' ,8X, 'EL' ,8X, 'TF' , 

00264:  &  3X,/, (3X,1P6E10.2) ,2X) 

00265:1109  CONTINUE 

00266:C 

00267:C 

00268:C  THETA  and  AMP  control  implicitness  in  time  stepping  and  iteration 
00269:C 

00270:  READ(5,*)  TIMLIM,  AMP,  THETA 
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00271: 
00272: 
00273: 
00274: 
00275: 
00276:' 
00277: 
00278: 
00279: 
00280: 
00281: 
00282: 
00283: 
00284: 
00285: 
00286: 
00287: 
00288: 
00289: 
00290: 
00291: 
00292: 
00293: 
00294; 
00295: 
00296: 
00297: 
00298: 
00299: 
00300: 
00301: 
00302: 
00303: 
00304: 
00305: 
00306: 
00307: 
00308: 
00309: 
00310: 
00311 : 
00312: 
00313: 
00314: 
00315: 
00316: 
00317: 
00318; 
00319: 
00320: 
00321 : 
00322:' 
00323:' 
00324:' 


WRITE(IO,706)  TIMLIM,  AMP,  THETA 
706  FORMAT (////, 15X Solution  Control  Parameters 21X , 

&  'TIMLIM' ,9X, 'AMP' ,7X, 'THETA' ,/,15X,lP3E12.3) 

Determine  ICS 

READ(5,*)  KIC 
IF(KIC.EQ.l )  GO  TO  60 
READ{5,*)  TEMPIC 

Initialize  TOLD 

DO  50  ND  =  1,NBND 

TOLD(ND)  =  TEMPIC 
50  CONTINUE 

GO  TO  99 

60  DO  70  ND  =  1,NBND 

READ(5,*)  NODE,  TEMPIC 
TOLD(NODE)  =  TEMPIC 
70  CONTINUE 

99  CONTINUE 

WRITE( 10, 733)  (K,TOLD(K) ,K=1,NBND) 

733  FORMAT( /////, 15X ,' Ini tial  Temperatures  at  Each  Node',//, 
&  3(I10,1P1E12.3)  ) 


Set  BC'S  and  time  step  values 
READ(5,*)  KBCT 

For  KBCT  =  1,  be  values  are  constant  over  time,  with 
KDT  =  1  for  uniform  time  steps,  or 
KDT  =  2  for  uniform  steps  in  SQRT(TIME) 

for  KBCT  =  2  be  values  over  time  and  time  steps  are  input 

IF(KBCT.EQ.2)  GO  TO  360 

Constant  EC's,  initialize  the  RHS  of  the  nonlinear  system  at  bc's 

DO  301  ND  =  1,NMBC1 

READ(5,*)  NODE,  TBCVAL 
R(NODE)  =  TBCVAL 
301  CONTINUE 

READ(5,*)  KDT 
IF(KDT.EQ.2)  GO  TO  310 

Uniform  time  steps,  calculate  DT(nt3) 


J*  ^ ^  ^ 


00325 

00326 

00327 

00328 

00329 

00330 

00331 

00332 

00333 

00334 

00335 

00336 

00337 

00338 

00339 

00340 

00341 

00342 

00343 

00344 

00345 

00346 

00347 

00348 

00349 

00350 

00351 

00352 

00353 

00354 

00355 

00356 

00357 

00358 

00359 

00360 

00361 

00362 

00363 

00364 

00365 

00366 

00367 

00368 

00369 

00370 

00371 

00372 

00373 

00374 

00375 

00376 

00377 

00378 


DTCNST  =  TIMLIM/NTS 
DO  305  IDT  =  1,NTS 
305  DT(IDT)  =  DTCNST 
GO  TO  399 

Uniform  steps  in  SQRT{TIME)(  calculate  DT(nts) 

310  DTSQT  =  SQRT( TIMLIM)/NTS 
DT(1)  =  DTSQT**2 
TIM2  =  DT(1) 

DO  320  IDT  =  2, NTS 
TIMl  =  TIM2 
TSQT  =  SQRT(TIM2) 

TNWSQT  =  TSQT  +  DTSQT 
TIM2  =  TNWSQT**2 
DT(IDT)  =  TIM2  -  TIMl 
320  CONTINUE 
GO  TO  399 


Read  in  be  values  and  associated  times;  with  number  of  intervening 
time  steps.  Set  flags  and  interpolate  between  nodes. 

360  CONTINUE 
IT  =  0 

TIMOLD  =  0.0 

DO  390  IBV  =  1,10000 

READ(5,*)  TIMBD,  NUMCHG,  NTSBD 
DTIM  =  (T1MBD-TIM0LD)/NTSBD 
TIMOLD  =  TIMBD 
DO  365  ND  =  1,NMBC1 

365  NBCFLG(ND)  =  -1 

N1  =  1 

DO  380  ND  =  l/NMBCl 

IF(ND.GT. NUMCHG)  GO  TO  366 
READ (5,*)  NODE,  TBCVAL 
TBCVLP  =  TOLD(NODE) 

NM  =  NGLNBC(NODE) 

IF(IT.NE.O)  TBCVLP  =  TBC(IT,NM) 

DTMP  =  (TBCVAL-TBCVLP) /NTSBD 
GO  TO  377 

366  CONTINUE 

DO  370  NM  =  N1,NMBC1 

IF(NBCFLG(NM).LT.O)  GO  TO  371 

370  CONTINUE 

371  N1  =  NM+1 

NODE  =  NBCNGL(NM) 

IF( IT.EQ.O)  GO  TO  374 
TBCVLP  =  TBC(IT,NM) 

ITM  =  IT-1 

IF(IT.EQ.l)  TBCVPP  =  TOLD(NODE) 

IF(IT.GT.l)  TBCVPP  =  TBC{ITM,NM) 


00379 

DTMP  =  (TBCVLP-TBCVPP)  •  DT  I  DT  (  I T  ) 

'  V 
^  i 

00380 

GO  TO  377 

V 

W 

00381 

374  TBCVLP  =  TOLD(NODE) 

00382 

DTMP  =  0.0 

1 

00383 

377  NBCFLG(NM)  =  1 

00384 

DO  378  ITB  =  l.NTSBD 

00385 

ITC  =  IT  +  ITB 

00386 

DT( ITC)  =  DTIM 

00387 

TBC(ITC,NM)  =  TBCi/LP  +  ITB*DTMP 

00388 

378  CONTINUE 

a 

00389 

380  CONTINUE 

h 

00390 

C 

00391 

IT  =  IT  NTSBD 

00392 

IF( IT.GE.NTS)  GO  TO  399 

00393 

390  CONTINUE 

00394 

C 

•C‘i 

00395 

399  CONTINUE 

K 

00396 

C 

00397 

C 

00398 

IF(KBVPRT.EQ.2)  GO  TO  440 

00399 

IF(KBVPRT.NE.l)  GO  TO  499 

00400 

WRITE( 10,741 ) 

00401 

741  FORMAT( //// , 15X , ' Time  and  DT  values  at  Each  Time  Step',//, 

00402 

&  11X,'TIME  STEP'  ,8X,  'TIME'  ,10X,  'DT'  ) 

I 

*» 

00403 

TIMCK  =  0.0 

Si 

00404 

DO  401  IT  =  1 , NTS 

00405 

TIMCK  =  TIMCK  +  DT(IT) 

00406 

WRITE( 10, 742 )  IT  ,  TIMCK , DT ( IT ) 

00407 

742  FORMAT! 120, 1P2E12. 3) 

00408 

401  CONTINUE 

1 

00409 

GO  TO  499 

00410 

C 

00411 

440  CONTINUE 

00412 

WRITE! 10,1743 ) 

00413 

1743  format! /////, 1 5X ,' Time  and  DT  at  Each  Time  Step,  with  ',/, 

00414 

&  20X ,' Boundary  Values  at  Each  Node') 

,s* 

00415 

TIMCK  =  0.0 

00416 

DO  460  IT  =  1,NTS 

00417 

TIMCK  =  TIMCK  +  DT!IT) 

004ia 

WRITE! 10, 743 )  IT , TIMCK , DT ! IT ) 

00419 

743  FORMAT!////, lOX, 'TIME  STEP  = '  , 14 , 5X ,  ' TIME  = '  , 1 P 1  El  2 . 3 , 5 X 

f 

00420 

4  'DT  =' ,1P1E12.3) 

' « 

00421 

WRITE! 10,745) 

00422 

745  FORMAT! /, 20X ,' Nodes  and  their  temperatures') 

1. 

00423 

IF!KBCT.EQ.2) 

00424 

&  WRITE!I0,744)  !NBCNGL! J) ,TBC! it, J ) , J=1 ,NMBC1 ) 

00425 

744  FORMAT! 3! 110 , 1P1E12.3) ) 

00426 

IF!KBCT.EQ.l) 

00427 

&  WRITE! 10, 744)  ! NBCNGL ! J ) , R ! NBCNGL ! J ) ) , J=1 , NMBCl ) 

*  /j 

00428 

460  CONTINUE 

00429 

C 

00430 

499  CONTINUE 

00431 

C 

00432 

C 

'it 


Input  convective  BC  stuff 


00433;C 
00434:C 

00435:  IF(KBC2.Eg.O)  GO  TO  44 

00436:  READ(5,*)  HH ,  TA 

00437:  WRITE( 10, 7744 )  HH,  TA 

00438:  7744  FORMAT (//, 20X ,' HH  = '  , 1  PI  El  2 . 3 , 4X , ' TA  =',1P1E12.3) 

00439;  44  CONTINUE 

00440:C 

00441:  IF(KPAUSE.EQ.O)  GO  TO  998 

00442:  WRITE( 1 , 7998 ) 

00443:  7998  FORMAT( ///////4X , ' OK  to  proceed  (  1  =  VES)') 

00444:  READd,*)  KOK 

00445:  WRITE( 1 , 8000 ) 

00446:8000  FORMAT( //////, 20X ,' Computing  Finite  Element  Solution...',/////) 
00447:  IF(KOK.EQ.l)  GO  TO  998 

00448:  WRITE( 1 , 7999 ) 

00449:  7999  FORMAT( //////////, lOX ,' User  mandated  STOP!') 

00450:C 

00451 :C  Close  all  files  and  stop  if  user  has  so  commanded 

00452:C 

00453:  CLOSE(5) 

00454:  CLOSE(IO) 

00455:  CLOSE(33) 

004  56  :C 

00457:  STOP 

00458:C 

00459:  998  CONTINUE 

00460:C 

00461  :C 

00462;C  Initialize  time  loop  parameters 
00463  :C 

00464;  TIME  =  0.0 

00465:  THETAM  =  THETA  -  1.0 

00466:  AMPM  =  1.0  -  AMP 

00467:  INEW  =  0 

00468:  ITRLMP  =  ITRLM+1 

00469:  KTPRNT  =  0 

00470:  KTPLOT  =  -1 

00471:C 

00472:C  Adjust  dt  and  tbc 
00473:C 

00474:  NTSP  =  NTS  +  1 

00475:  DT(NTSP)  =  DT(NTS) 

00476:  DO  999  NM=1,NMBC1 

00477:  999  TBC(NTSP,NM)  =  TBC(NTS,NM) 

00478:C 

00479:  IF( I2DPLT.NE.0)  WR ITE ( 33 , 7099  )  IZDPLT 

00480:  IFdZDPLT.EQ.O)  IZDPLT  =  -1 

00481:  7099  FORMAT(I5) 

00482:  C 

00483:C  Initialize  ANS,  TLAST 

00484  :C 

00485:  DO  989  ND  =  1,NBND 

00486:  ANS(ND)  =  TOLD(ND) 


00487:  TLAST(ND)  =  TOLD(ND) 

00488:  989  CONTINUE 

00489:C 

00490:  WRITE( 10, 7001 ) 

00491:  7001  FORMAT ( 5 (/////), 20X Ca Icuia ted  Temperature  Solution',/////) 

00492:C 

00493:C 

00494:C 

00495:C 

00496:C  Time  Loop 

00497:C 

00498:C 

00499:  DO  1000  IT  =  1,NTSP 

00500:C 

00501:  TIMOL  =  TIME 

00502:  TIME  =  TIME  +  DT(IT) 

00503:  ITER  =  0 

00504:  KTPRNT  =  KTPRNT  +  1 

00505:  KTPLOT  =  KTPLOT  +  1 

00506:C 

00507:  IF(KBCT.EQ.l )  GO  TO  930 

00508:  DO  920  NM  =  1,NMBC1 

00509:  NODE  =  NBCNGL{NM) 

00510:  R(NODE)  =  TBC(IT,NM) 

00511:  920  CONTINUE 

00512:  930  CONTINUE 

00513:C 

00514:C  Call  to  routine  which  performs  element  by  element  integrations  and 

00515:C  returns  coefficient  matrix  A  and  r.h.s  vector  R 

00516:C 

00517:  901  CALL  TR I  NT { A , NBND , NBW , NDBC , NODEL , NELS , XG , VG , 10 , NDG , DT ( IT )  , R , 

00518:  &  CNDF,CNDU,HETF,HETU,TF, THETA, THETAM,T, TOLD, TLAST, EL, ITER, 

00519:  &  ITRLMP,TIMOL, AMP, AMPM, IZDPLT , KTPLOT , ANS , IT , MEDIA , MTYPE , HH , TA ) 

00520:  C 

00521:C  Solve  algebraic  system  of  equns 
00522:C 

00523:  CALL  BANSAL ( A , NBW , NBND , NBND, R , ANS , I  NEW , KERR , 10 ) 

00524:C 

00525:  IF( ITER. EQ. ITRLM)  GO  TO  950 

00526:  ITER  =  ITER+1 

00527:  GO  TO  901 

00528:C 

00529:  950  CONTINUE 

00530:C 

00531:  IF(KTPRNT.NE.KPRTVL)  GO  TO  990 

00532:  KTPRNT  =  0 

00533:C 

00534:C  Print  results 

00535:C 

00536:  WRITE( 10, 716  )  TIME,  IT 

00537:  716  FORMAT{ // ,  25X  ,  'TIME  = '  , 1  Pi E 1 2 . 3 , 4X  ,  ' T IME  STEP  =',I6) 

00538:  WRITE! 10,713)  (  K  ,  ANS ( K ), K= 1 , NBND  ) 

00539:  713  FORMAT! // , 3!6X , ' NODE’  ,8X,  'TEMP'  ),/,3!5X,I5,lPlEl2.3)) 

00540:  990  CONTINUE 


00541  :C 

00542:C  Update  values  of  TLAST  and  TOLD 
00543:C 

00544:  DO  800  ND  =  1,NBND 

00545:  TLAST(ND)  =  TOLD(ND) 

00546:  800  TOLD(ND)  =  ANS(ND) 

00547  :C 
00548:C 

00549:  1000  CONTINUE 
00550:C 

00551:  FIL  =  1.0 

00552:  VRYNEG  =  -1.1E30 

00553:  IF( IZDPLT.GT.O)  WR ITE ( 33 , 7007 )  VRYNEG, FIL , FIL , FIL 

00554:  7007  FORMAT ( 1P4E1 1 . 2 ) 

00555:C 

00556:C 

00557:  1111  CONTINUE 
00558:C 

00559:C  Close  all  files  and  stop 

00560:C 

00561:C 

00562:  CLOSE(5) 

00563:  CLOSE(IO) 

00564:  CLOSE(33) 

00565:C 

00566:C 

00567:  RETURN 

00568;  END 

00569:C 

00570:C 

00571:C 

00572:C 

00573:C 

00574:C 

00575:  SUBROUTINE  TRINT ( A , NBND, NBW , NDBC , NODEL , NELS , XG , YG , 10, NDG , DT , 

00576 :  &  R , CNDF , CNDU , HETF , HETU , TF , THETA , THETAM , T , TOLD , TLAST , EL , ITER , 

00577:  &  ITRLMP,TIMOL, AMP, AMPM, I2DPLT,KTPLOT,ANS, IT,MEDIA,MTYPE, HH,TA) 

00578:C 

00579;C  Subroutine  to  perform  finite  element  integrations  and  assemble  banded 
00580:C  global  matrix  for  linear  triangles 

00581:C 
00582  :C 

00583:  COMMON  /SMALL/  ZERO 

00584:C 

00585;  DIMENSION  A{NBND, NBW) ,  NDBC(NBND),  NODEL ( NELS , 3 ) ,  XG(NBND) 

00586:  DIMENSION  YG(NBND),  X(3),  Y(3),  R(NBND),  T(NBND),  TOLD(NBND) 

00587;  DIMENSION  TLAST(NBND),  ELA(3),  ELB(3),  TE(3),  ANS(NBND),  EL(IO) 

00588:  DIMENSION  TF(IO),  HETU ( 10 ) , HETF ( 10 ) 

00589:  DIMENSION  CNDF ( 10 ), CNDU ( 10 ), MTYPE ( NELS ) 

00590:C 

00591:C  Identifies  a  first  pass 
00592  :C 
00593: 

00594: C 


IF( IT* ( ITER+1 )  .NE,  1)  GO  TO  101 


00595:C  The  next  computations  eliminate  the  problems  that  could 
00596;C  start  if  the  maximum  initial  temperature  was  zero 
00597  :C 

00598:  TMAX  =  TOLD(l) 

00599:  TMIN  =  TOLD{l) 

00600:  DO  100  ND  =  2 , NBND 

00601:  TM  =  TOLD(ND) 

00602:  IF(TM  .GT.TMAX)  TMAX  =  TM 

00603:  IF(TM  .LT.TMIN)  TMIN  =  TM 

00604:  100  CONTINUE 

00605:  DELTF  =  0.0025*(  TMAX  -  TMIN  ) 

00606:  TMAX  =  ABS(TMAX) 

00607:  IF(ABS(TMIN) .GT.TMAX)  TMAX  =  ABS(TMIN) 

00608;  IF(DELTF.LT.  0.0025*TMAX)  DELTF  =  0.0025*TMAX 

00609:  101  CONTINUE 

00610:C 

00611:C  Update  the  T(ND)  used  to  calculate  the  matrix  A 
00612:C 

00613:  DO  105  ND  =  1,NBND 

00614:  T(ND)  =  AMP*ANS{ND)  +  AMPM*TOLD ( ND ) 

00615:  105  CONTINUE 

00616:C 

00617:  DO  10  ND  =  1,NBND 

00618:C 

00619:C  Initialize  the  matrix  A,  the  r.h.s.  vector  R 
00620:C 

00621:  DO  5  J  =  1,NBW 

00622;  5  A(ND,J)  =  0.0 

00623:  1F(NDBC(ND) .NE.l )  GO  TO  9 

00624:  A(ND,NDG)  =  1.0 

00625:  GO  TO  10 

00626:  9  R{ND)  =  0.0 

00627:  10  CONTINUE 

00628:C 

00629:C 

00630:C  ELEMENT  LOOP.  This  is  the  one  that  does  the  quadrature/ 

00631:C  and  assembles  A  and  R 

00632:C 

00633:  DO  1000  NLMT  =  1,NELS 

00634;C 

00635:C 

00636:C  Initialize  flags 

00637:C  KPHS  =  -3  :  Element  frozen 

00638:C  KPHS  =  +3  :  Element  unfrozen 

00639;C  KFRZ  =  1  :  Phase  change  element 

00640:  KFRZ  =  3 

00641:  KPHS  =  0 

00642;  RF  =  1.0 

00643:C 

00644:C  Identify  material  type  of  this  element 
00645  :C 

00646:  II  =  MTYPE(NLMT) 

00647:C 

00648:C  The  next  loop  works  with  each  element  to  determine  if  it  is 
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a  phase  change  element  or  not  and  then  does  necessary 
computations  to  enable  calculation  of  phase  change  effect 


DO  400  K  =  1,3 

NK  =  N0DEL( NLMT, K) 

X(K)  =  XG(NK) 

Y(K)  =  YG{NK) 

TE(K)  =  T(NK) 

This  section  accounts  for  the  possibility  that  the  node 

temperature  equals  the  phase  change  temperature.  If 

it  does,  it  is  shifted  in  the  direction  of  its  last  recorded  value 

TFM2  =  TF(Il)  -  DELTF 
TFP2  =  TF( II )  +  DELTF 

IF  ( (TE( K) .GT.TFMZ ) .AND. (TE( K ) . LT.TFPZ ) )  THEN 
IF  (TLAST(NK  )  .GT.TF( II )  )  THEN 
TE(K)  =  TF(I1)+DELTF 
ELSE 

TE(K)  =  TF(Il)  -  DELTF 
ENDIF 
ENDIF 

Record  frozen  or  unfrozen  status  of  node 

IF(TE(K) .LT.TF( II ) )  KPHS  =  KPHS  -  1 
IF(TE(K ) .GT.TF( II ) )  KPHS  =  KPHS  +  1 

Identify  next  node  around  the  element 

KP  =  K+1 

IF(KP.EQ.4)  KP  =  1 
NP  =  NODEL(NLMT,KP) 

Bypass  the  rest  of  the  loop  if  no  phase  change 

IF( ( (T(NP) .GT.TF( II ) ) .AND. ( T ( NK ) . GT . TF ( I 1 ) ) ) . OR . 

&  ( (T(NP) .LT.TF( II ) ) .AND. (T(NK ) . LT.TF( II ) ) ) )  GO  TO  400 

IF(KFRZ.EQ.3)  NPL  =  NP 
IF(KFRZ.EQ.3)  GO  TO  380 
NDF  =  NPL  -  NK 
IF(NDF.NE.O)  GO  TO  360 

(XC,YC)  is  vertex  of  triangular  sub-element 

NDl  =  NK 
XC  =  X(K) 

YC  =  Y(K) 

GO  TO  380 
360  NDl  =  NP 

XC  =  X(KP) 

YC  =  Y(KP) 
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380  CONTINUE 

C 

KFRZ  =  KFRZ-1 
C 

C  The  next  computations  are  for  the  phase  change  element/ 

C  computing  the  coordinates  of  the  phase  change  isotherm/ 

C  and  the  values  of  the  basis/weight  function  at  the  isotherm 
C  ( ELA  and  ELB ) . 

C 

C  The  next  few  lines  interpolate  to  find  the  phase 
C  change  isotherm 
C 

DLTF  =  ABS( (TF( Il)-T(NK) )/(T(NP)-T(NK) ) ) 

DELX  =  DLTF* (XG(NP)-XG{NK)  ) 

DELY  =  DLTF*( YG(NP)-YG(NK) ) 

IF(KFRZ.EQ.l  )  GO  TO  390 
XB  =  XG(NK)  +  DELX 
YB  =  YG(NK)  +  DELY 
ELB(KP)  =  DLTF 
ELB(K)  =  1.0  -  DLTF 
KPP  =  KP+1 

IF(KPP.EQ.4)  KPP  =  1 
ELB(KPP)  =  0.0 
GO  TO  400 

390  XA  =  XG(NK)  +  DELX 

YA  =  YG(NK)  +  DELY 
EL''(KP)  =  DLTF 
ELA(K)  =  1.0  -  DLTF 
KPP  =  KP  +  1 
IF{KPP.EQ.4)  KPP  =  1 
ELA(KPP)  =  0.0 
400  CONTINUE 

C  Update  TLAST 

DO  102  ND  =  1/NBND 
102  TLAST(ND)  =  T(ND) 

IF(KPHS.NE.3)  GO  TO  410 
C 

C  The  next  assignments  are  for  an  unfrozen  element 
C 

HET  =  HETU(Il) 

COND  =  CNDU(Il) 

GO  TO  499 

410  IF(KPHS.NE.-3)  GO  TO  420 

C 

C  The  next  assignments  are  for  a  frozen  element 
C 

HET  =  HETF(Il) 

COND  =  CNDF(Il) 

GO  TO  499 
C 

C  Calc  area  of  triangular  sub-element 


A1  =  XB*YC  -  XC*YB 
A1  =  ABS(Al/2.0) 


XA*YC  +  XC*YA  +  XA*YB  -XB*YA 


S'  .v  .s  s  .s 
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00757;C 

00758:C  Assign  frozen  and  unfrozen  conduct! vi  ties 

00759:C  to  appropriate  parts  of  element 

00760  :C 

00761:  COND  =  0.0 

00762:  VKl  =  CNDU(Il) 

00763:  VK2  =  CNDF(Il) 

00764:  IF{T(ND1 ) .GT.TF( II ) )  GO  TO  498 

00765:  VKl  =  CNDF{I1) 

00766:  VK2  =  CNDU(Il) 

00767:  498  CONTINUE 

00768:  499  CONTINUE 

00769:C 
00770:C 

00771:C  Calculate  the  area  of  the  element  using  cross  product 
00772:C 

00773:  AE  =  X{2)*Y(3)  -  X(3)*Y(2)  -  X(1)*Y(3) 

00774:  &  +  X(3)*Y(1)  +  X(1)*Y(2)  -  X(2)*Y(1) 

00775:  AE  =  AE/2.0 

00776:  C 

00777:C 

00778:C  Assign  conductivity  factor  for  non-phase  change  element 

00779:C 

00780:  C4A  =  0.25*COND/AE 

00781:C 

00782:  IF{KFRZ.EQ.l)  GO  TO  450 

00783:  BCF  =  HET  *  AE/(12.0*DT) 

00784:  GO  TO  4450 

00785; C 

00786:C 

00787:C  This  next  assignment  accounts  for  the  different 
00788:C  sensible  heat  in  the  unfrozen  and  frozen  portion 
00789:C  of  the  phase  change  element. 

00790:C 

00791:  450  CONTINUE 

00792:  IF  (T(NDl) .GT.TF(Il) )  THEN 

00793:  Cl  =  HETU(Il) 

00794:  C2  =  HETF(Il) 

00795:  ELSE 

00796:  Cl  =  HETF(Il) 

00797:  C2  =  HETU(Il) 

00798:  END  IF 

00799:  BCF  =  (Cl*Al  +  C2*(AE  -  A1 ) ) / ( 1 2 . 0*DT ) 

00800:C 

00801:C  Calculations  on  element  with  phase  change  to  account  for  the 

00802:C  latent  heat  effect  discontinuity 

00803:C 

00804:C  Calculate  length  of  tau 

00805:  TAU  =  SQRT(  {XB-XA)**2  +  (YB-YA)**2  ) 

00806:C  The  magnitude  of  the  gradient  of  T  is  computed  using  the  x  and  y 
00807:C  derivatives  of  the  equation  of  the  plane  T(x;y)  =  a  +  bx  +  cy / 

00808:C  i.e./  expressing  a  and  b  in  terms  of  x(i)»  y(i)»  &  T{i)  using  Cramer's 
00809:C  rule.  The  2*AE  is  the  determinant  used  in  Cramer's  rule. 

00810:  SUMDX  =  TE(l)*(Y(2)-y(3) )  +  TE( 2 ) * ( Y( 3 ) -Y{ 1 ) ) 
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&  +  TE(3)*(Y(1)-Y(2) ) 

SUMDY  =  TE(3)*(X{2)-X{1) )  +  TE ( 2 ) * ( X ( 1 ) -X ( 3 ) ) 

&  +  TE(1)*(X(3)-X(2)) 

DTDS  =  SQRT(  SUMDX**2  +  SUMDY** 2 ) / ( 2 .0*AE ) 
ELDSDT  =  EL( II )/{DTDS*DT) 

C 

C  Store  2 . d. isotherm  locations  for  plotting 
C 

IF(KTPLOT.EQ- IZDPLT.AND.ITER.EQ.O)  THEN 
WRITE(33/7001)  XA,YA,XB,YB 
7001  F0RMAT(1P4E13.4) 

ENDIF 

C 

4450  CONTINUE 

C 

C  ELEMENT  NODE  ROW  LOOP 
C 

DO  600  1=1,3 
C 

IG  =  NODEL( NLMT, I ) 

IF(NDBC( IG) .EQ.l )  GO  TO  600 
C 

C  Treat  convective  BC's 

C 

IF(NDBC( IG)  .NE.  2 )  GO  TO  599 
DO  510  J  =  1,3 
IF(J.EQ.I)  GO  TO  510 
JG  =  NODEL(NLMT,J) 

IF(NDBC(JG) .NE.2)  GO  TO  510 
DDX  =  X{I)  -  X(J) 

DDY  =  Y(I)  -  Y(J) 

DD  =  SQRT(DDX**2  +  DDY**2) 

R(IG)  =  R(IG)  +  HH*TA*DD/2.0 
&  +  THETAM*HH*DD  *  ( TOLD ( IG ) /3 . 0  +  TOLD ( JG ) /6 . 0 ) 

A(IG,NDG)  =  A(IG,NDG)  +  THETA* HH* DD/ 3 . 0 
JGBD  =  JG  -  IG  +  NDG 

A(IG,JGBD)  =  A(IG,JGBD)  +  THETA* HH*DD/6 . 0 
510  CONTINUE 
599  CONTINUE 
C 

IJ  =  I+l 

IF{IJ.EQ.4)  IJ  =  1 
IK  =  IJ+1 

IF(IK.EQ.4)  IK  =  1 
C 

C  ELEMENT  NODE  COLUMN  LOOP 
C 

DO  500  J  =  1,3 
JI  =  J-1 

IF(JI.EQ.O)  JI  =  3 
JK  =  J  +  1 
IF(JK.EQ.4)  JK  =  1 
JG  =  NODEL(NLMT, J) 

JGBD  =  JG  -  IG  +  NDG 
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C  This  term  (DRV)  is  part  of  the  element  quadrature 

DRV  =  (Y(IJ)-Y(IK) )*(Y(JK)-Y(JI)  )  +  ( X ( IK  ) -X { I J  )  ) 

&  *(X(JI)-X(JK) ) 

BIJ  =  BCF 

IF(I.EQ.J)  BIJ  =  2.0*BCF 
IF(KFRZ.NE.l )  GO  TO  480 
C 

C  For  the  phase  change  elements 
C 

BIJ  =  BIJ  +  ELDSDT*TAU*  (  2 . 0* ( ELB ( I ) -ELA { I ) ) * ( ELB ( J ) -ELA ( J ) ) 
&  +  3.0*ELB{I)*ELA(J)  +  3 . 0*ELB { J  )  *ELA ( I )  )/6.0 

C 

C  The  next  term  distributes  the  appropriate  conductivity 
C  over  the  correct  proportions  in  the  phase  change  element. 

C 

CIJ  =  VK2*DRV/(4.0*AE)  +  ( VK1-VK2 ) *DRV* Al/ ( 4 . 0* AE* AE ) 

GO  TO  495 
C 

C  For  the  other  elements 
C 

480  CIJ  =  C4A*DRV 

C 

C  Load  the  matrix  A  and  the  vector  R  for  the  system  solver 
C 

C  The  THETA  parameter  controls  numerical  implicitness 

C 

495  A(IG,JGBD)  =  A(IG,JGBD)  +  BIJ  +  THETA*CIJ 

R(IG)  =  R(IG)  +  (BIJ  +  THETAM*CIJ)*  TOLD(JG) 

C 

500  CONTINUE 

600  CONTINUE 

C 

1000  CONTINUE 


IF(KTPLOT.NE. IZDPLT  .OR.  ITER.NE.O)  GO  TO  1010 
WRITE(IO,7070)  TIMOL 

7070  FORMAT(//,12X, 'Plotting  data  stored  for  TIME  =',1P1E12.3) 
KTPLOT  =  0 
VRYLGE  =  1.1E30 
VITM  =  IT-1.0 

WRITE(33,7701)  VRYLGE,  VITM,  TIMOL,  TIMOL 
7701  FORMAT( 1P4E11.2) 

1010  CONTINUE 


Determine  "ZERO",  i.e.  a  magnitude  below  which  significance  is  lost 

IF(IT*(ITER+1)  .NE.  1)  GO  TO  1101 

AMAX  -  0.0 

DO  1100  I  -  1,NBND 

DO  1100  J  -  1,NBW 

ABSVAL  -  ABS(A(I,J) ) 

IF(ABSVAL.GT.AMAX)  AMAX  «  ABSVAL 
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CONTINUE 

ZERO  =  AMAX  *  0.0001 


CONTINUE 

RETURN 


C 

C 

C 

C  The 


following  routine  solves  the  system  Ax=B 


SUBROUTINE  BANSAL  { A , NBW , NEQ , NMAX , B , X , INEW , KERR , 10 ) 


COMMON  /SMALL/  ZERO 

DIMENSION  A(NMAX,NBW) ,B(NMAX) ,X(NMAX) 

INITILIZE  COUNTERS 

KERR=0 

NBWT=NBW 

N=NEQ 

NB2=  NBWT/2  +1 
NB1=  NB2-1 
N1=N-1 
NP1=N+1 
NB1M1=  NBl-1 
DO  10  1=1, N 
X(I)=B(I) 

CONTINUE 

IF( INEW.NE.O)GO  TO  350 

FORWARD  ELIMINATION  OF  MATRIX 

K1=NB1 
IM1  =  1 

DO  300  1=2, N 

K1=MIN0(K1+1,N) 

K  =  NB1 

AII=A(IM1,NB2) 

IF( (All. LT.ZERO) .AND. (All. GT. -ZERO) )  GO  TO  1850 
DO  200  J=I,K1 

CX=  A(J,K)/  All 
Ll=NB2  +1 
KPB=  K  +NB1 
KP1=  K+1 

DO  100  L=KP1,KPB 

A(J,L)=  A(J,L)  -CX*  A(IM1,L1) 

L1=L1+1 


V.'.’Ls' 
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A(J,K)=  CX 
X( J)=X( J)-X( IMl )*CX 
K=  K-1 
200  CONTINUE 

300  IM1=I 

GO  TO  550 
C 

C  OPERATE  ON  THE  NEW  RHS  -IS  DONE  AFTER  THE  FIRST  PASS 

C 


350 


400 

500 

C 

C 

C 

550 


IN=NB1 
IM1  =  1 

DO  500  1  =  2, N 

IN=MIN0(IN+1/N) 

Jl=  NBl 
DO  400  J=I , IN 

X(J)=  X(J)  -X(IMl)*  A(J,K1) 
K1=K1-1 
CONTINUE 
IM1  =  I 

BACK  WARD  SUBSTITUTION 

X{N)=  X(N)/A(N,NB2) 

DO  700  1=1, N1 
SUM=0.0 
L=  NB2  +1 
K=  N-I 

Kl=  MIN0(N1,  K+NBIMI) 

DO  600  J=K,K1 

SUM=SUM  +  A(K,L)*  X(J+1) 

L  =  L  +  1 


600  CONTINUE 

X(K)  =  (  X(K)-SUM)/  A{K,NB2) 

700  CONTINUE 

RETURN 

1850  WRITE(IO,9000)  IMl 
KERR=1 
RETURN 

9000  FORMAT(1HO,120(1H*)/10X, 'ERROR  FROM  BANSAL- ' , 15 , 13H  DIAGONAL  ELE, 
1  20HMENT  REDUCED  TO  ZERO/ IX , 1 20 ( 1 H* )  ) 

END 
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